AD-A238  858 

llllllllll 


® 


Grant  AFOSR  86-0337 

Dire 

F!.rC ' 

J  U  1_  lit  o  i  s  J  i  L" 

BASIC  PROCESSES 

OF 

M  C  ** 

PLASMA  PROPULSION 

Herbert  0.  Schrade 
P.  Christian  Sleziona 
Helmut  L.  Kurtz 


A  :  'j#.i  iliic  t  ir 
*  '  r;."'.4l 

'  n.-» 

J  .sici.*  *.  c  >t  U4i. 


April  1991 


Final  Scientific  Report 
1  August  1986  -  31  July  1990 


•  lu/or 


Approved  for  public  release;  distribution  unlimited 
Prepared  for  the 

Air  Force  Office  of  Scientific  Research 
Bolling  AFB,  DC  20332-6448 
through  the 

European  Office  of  Aerospace  Research 
and  Development 
223/231  Old  Marylebone  Road 
London,  NW1  5th 
England 


91-05931 

lllllllll 


•1  7  23  012 


DISCLAIM!!  NOTICE 


THIS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTIC  CONTAINED 
A  SIGNIFICANT  NUMBER  OF 
COLOR  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY  ON  BLACK 
AND  WHITE  MICROFICHE. 


Preface 


The  work  descriped  in  this  final  report  has  been  sponsored  by  the  Air  Force  Office  of 
Scientific  Research  through  the  European  Office  of  Aerospace  Resarch  and 
Development  in  London  under  Grant  AFOSR  86  -  0337.  The  research  work  period 
herein  is  from  August  1,  1986  to  July  31,  1990. 


I2£R£E] 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No  0704  0188 


'a  REPORT  security  Classification 

ib  restrictive  markings 

2.  SECURITY  CLASSIFICATION  AUTHORITY 

.... 

.  2b  DECLASSIFICATION /DOWNGRADING  SCHEDULE 

'  i  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

S.  M0> STORING  ORGANIZATION  REPORT  NUMBER(S) 

6*  NAME  OF  PERFORMING  ORGANIZATION 

Inst.  f.  Raumfahrtsysteme 
University  of  Stuttgart 

6b  OFFICE  SYMBOL 

(if  spplicsble ) 

IRS 

7a.  NAME  OF  MONITORING  ORGANIZATION 

EOARD ,  London 

6c.  ADDRESS  (City,  State,  end  HP  Cod*) 
Pfaf fenwaldring  31 

7000  Stuttgart  80 


Germany 


7b  A0MEft(Cfty,  SUM.  anrfWCodt) 

tfOSWW  *7  Mb'*  MMMLAfltfl 

soiling  jot  #c .  wm  hot 


8a  NAME  OF  FUNDING  /SPONSORING 
RGAN1ZATION 


8b  OFFICE  SYMBOL 

(If  4pfilk»b/f) 


•oiling  JOT  0CJ  I0332-6448] 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

AFOSR  Grant  86-00337 


to  SOURCE  OF  FUNOING  NUMBERS 


PROGRAM 

LEMENTJ 


PROJECT 

NO 


TASK 

NO 


WORK  UNIT 
I  ACCESSION  NO 


1 1  TITLE  (Include  Security  Classification) 


12 


13a 


16 


Basic  Processes  of  Plasma  Propulsion  (not  classified) 


PERSONAL  AUTHOR(S) 

Herbert  0.  Schrade ,  p. Christian  Sleciona,  Helmut  L.  Kurtz 

TYPE  OF^ REPORT  13b  TIME  COVERED  14.  DATE  OF  REPORT  (Yesr.  Month.  Dey) 

Final  FRpM  Auq.86  TOjuly90  199 T,  Jan,|  15 

SUPPLEMENTARY  NOTATION 


[15  PAGE  COUNT 


17  COSATI  CODES  i 

18  SUBJECT  TERMS  (Cont/nu*  on  rtv.ru  if  neettury  tnd  idtntify  by  black  number) 

FIELD 

GROUP 

SUB-GROUP 

19  ABSTRACT 

(Continue  on 

reverse  if  necesssry 

tnd  identify  by  block  number ) 

In  order  to  improve  perfomance,  lifetime  and  reliability  of  plasma  thrusters  one  has  to  understand  and 
asses  the  fundamental  processes  and  problem  areas  like  nonequilibrium  magneto  plasma  dynamics  flows 
in  electric  discharges,  their  stability  and  electrode  erosion  effects.  In  order  to  asses  the  flow,  pressure 
and  density  or  temperature  fields  in  plasma  thruster  devices  several  theoretical  model  calculations  have 
been  developed  and  applied  to  different  thruster  configurations.  The  results  of  these  calculations  are 
compared  with  corresponding  experimental  arc  devices  and  both  agree  well  with  each  other.  Concerning 
the  plasma  stability  work,  a  new  explanation  of  the  so  called  "Onset  Phenomenon"  is  presented.  It 
is  based  on  a  run-away  Joule  heating  effect  caused  by  the  self  magnetic  field  of  a  current  carrying 
plasma  channel  which  drastically  decreases  the  radial  heatconduction  losses  due  to  the  electrons.  The 
dependency  of  the  electron  heatflux  vector  and  heatconduction  coefficient  in  a  transverse  magnetic  field 


SECURITY  C.ASaFtCA^ON  OF  THt$  PAGE 


Previous  editions  ere 


k'sen- 


19  ABSTRACT  (Continue  on  reverse  if  netetssry  snd  identify  by  block  number) 

lias  been  derived  by  means  of  a  perturbation  approach  of  the  Maxwell  distribution.  Cathode  erosion 
measurements  have  been  conducted  by  means  of  a  fully  automatically  working  test  rig.  The  cathode 
sample  is  repeatably  charged  by  a  fairly  rectangular  current  pulse  of  1400  A  which  lasts  about  2  ms. 
Erosion  rate  measurements  on  thoriated  ( 2%ThO 2)  cathode  samples  are  reported  on  and  explained  by 
a  refined  cathode  spot  theory. 


Contents 


1  Introduction  2 

2  MPD  -  Code  Development  4 

2.1  The  Nozzle  Type  MPD  Thruster  .  4 

2.1.1  Flow  field  code .  5 

2.1.2  Discharge  code .  6 

2.1.3  Electron  temperature  code .  7 

•  2.1.4  Solution  of  the  Equations  .  8 

2.1.5  Discussion  of  Results .  9 

2.1.6  Conclusions .  11 

2.2  The  Cylindrical  MPD  Thruster .  12 

2.2.1  Solution  of  the  Equations  .  16 

2.2.2  Numerical  Results  .  16 

2.2.3  Comparison  with  experimental  results  .  30 

3  Plasmastability  34 

4  Electrode  Effects  45 

4.1  Cathode  Experiments .  45 

4.2  Measurements .  45 

4.3  Conclusion  and  Explanation  (qualitative) .  53 

1 


5  List  of  Reports  and  Publications  57 

6  Acknowledgements  60 

A  Appendices  62 

A.l  Influence  of  the  Magnetic  Field  on  the  Heatflux  within  an  Electrically 

Charged  Plasmacomponent .  63 

A.2  Determination  of  the  Temperatureprofile  within  the  Constrictor .  70 

Bibliography  78 


II 


List  of  Figures 

1.1  Different  MPD-Thruster  Devices  a)  The  Nozzle  Type  Thruster  b)  The 

Cylindrical  Type  Thruster .  3 

2.1  Thruster  with  the  Calculation  Areas.  All  length  Coordinates  are  in  [mm]  18 

2.2  Magnetic  Reynold’s  Number  Distribution  of  the  Three  Different  Throat 

Geometries  for  2  kA  and  a  Mass  Flow  of  0.8  g/s  Argon.  As  Characteristic 
Length  the  Nozzle  Diameter  was  Chosen .  19 

2.3  Calculated  Current  Contour  Lines  of  the  Three  Different  Throat  Geome¬ 
tries  for  2  kA  and  a  Mass  Flow  of  0.8  g/s  Argon .  20 

2.4  Electron  Temperature  Distribution  of  the  Three  Different  Throat  Geome¬ 
tries  for  2  kA  and  a  Mass  Flow  of  0.8  g/s .  21 

2.5  Calculated  and  Measured  Electron  Temperature  Profile  for  the  Thruster 
with  24  mm  Throat  Diameter  at  a  Current  Level  of  2  kA  and  a  Mass  Flow 

of  0.8  g/s .  22 

2.6  Potential  Lines  of  the  Three  Different  Throat  Geometries  for  2  kA  and  a 

Mass  Flow  of  0.8  g/s .  23 

2.7  Calculated  and  Measured  Discharge  Voltage .  24 

2.8  Comparison  of  Calculated  and  Measured  Thrust  for  a  Mass  Flow  of  0,8 

g/s  Argon  ct  =  1 .  25 

2.9  Calculated  and  Measured  Arc  Chamber  Pressure .  26 

2.10  Density  Maps  a)  with  and  b)  without  the  Pinch  Effect  of  the  DT2  with 

Argon  Discharge  at  2  kA  and  a  mass  Flow  of  0.8g/s .  26 

2.11  Calculation  Grid  of  the  ZT2-IRS  Thruster .  27 

2.12  Computation  Scheme .  27 

2.13  Magnetic  Reynold’s  Number  Distribution .  28 


III 


2.14  Calculated  Current  Contour  Lines .  28 

2.15  Measured  a)  and  Calculated  b)  Current  Contour  Lines .  29 

2.16  a)  Heavy  Particles  and  b)  Electron  Temperature  Contours .  32 

2.17  Density  Distribution .  33 

2.18  Velocity  Vector  Distribution .  33 

3.1  Illustration  of  the  Self  Magnetic  Field  within  a  Plasma  Thruster .  35 

4.1  Scheme  of  Cathode  Testrig .  46 

4.2  Testrig  for  Cathode  Experiment .  47 

4.3  Typical  Current  vers.  Time  Curve .  48 

4.4  Measured  Erosionrates  in  as  Funktion  of  the  Accumulated  Total  Elec¬ 

tric  Charge  in  [As]  for  Different  Ambient  Pressures.  Each  Measuring  Point 
is  the  Result  of  50  Current  Pulses. (Cathode:  Thoriated  Tungsten  at  Room 
Temperature) .  49 

4.5  Erosion  Rate  as  Funktion  of  the  Accumulated  Charge  (  Number  of  Current 

Pulses) .  50 

4.6  Asymptotic  Erosion  Rate  as  Funktion  of  the  Ambient  Pressure  of  a  Tungs¬ 
ten  Cathode  in  a  Nitrogen  Atmosphere  Charged  by  a  Rectangular  Pulses 

of  1400  A  for  2  ms  (Pulsfrequency  =  1-2  min~l) .  50 

4.7  Massloss  Due  to  Pure  Heating  in  a  Nitrogen  and  Argon  Atmosphere  (p  < 

100Pa) .  51 

4.8  Temperature  of  the  heated  probe  (e  =  0.45)  as  Funktion  of  the  Heating 

Power .  52 

4.9  Scheme  of  a  Spot  Discharge .  53 

4.10  Illustration  of  the  Plasmaflow  around  the  Crater  Rim  with  Prandtl- Meyer 

Deflection  Angle, 6 .  54 

A.l  Illustration  of  the  Particle  Motion  within  a  Magnetic  Induction  Field,  B, 
normal  to  the  Gradient,  V /,  :  (a)  towards  V /,;  (b)  V,±  opposite  to 
Vfi  (B  normal  to  plane) .  66 

A.2  Temperature  Profile  for  Different  Grades  n  of  the  Current  Density  Para- 

boloides  within  the  Arc  Core .  76 


IV 


A. 3  Plot  of  7  as  Funktion  of  Tc  and  p  as  Parameter  of  Hydrogen  as  Propellant  77 


A.4  Maximum  Electrontemperature  within  an  Hydrogen  Arc  Core  as  Function 
ot  the  Quotient  -jr  (arc  current  divided  by  the  core  radius)  for  differnt 
pressure  p,  n  =  4,  £  =  0 . 


77 


V 


List  of  Tables 


2.1  Comparison  of  Numerical  and  Experimental  Results  for  2  g/s  Argon  at 

6000  A .  30 

3.1  J2  according  to  eq. (3. 12)  and  (3.13) .  37 

3.2  0  for  Hydrogen .  38 

3.3  0  for  Argon .  38 

3.4  0’  for  Hydrogen .  41 

3.5  0‘  for  Argon .  42 

3.6  w,, r,. -Values  calculated  for  different  Onset  conditions  at  n=2,  4, 10  for  Argon  44 

3.7  w*Ts-values  calculated  for  different  Onset  conditions  at  n=2,  4,  10  for 

Hydrogen  .  44 

4.1  Erosionrate  caused  by  arcing  only  at  the  surface  of  a  heated  tungsten  (2% 
thoriated)  cathode  in  a  AV  atmosphere  (<  100Pa)  surface  temperature: 

«  1500 K .  52 

4.2  Prandel-  Mayer  Deflection  f  for  various  pressure  ratios  jja  and  k-  values  .  56 

A.  I  J0{n,  4)  .  72 

A. 2  J\(n,  O  .  73 


1 


1 


Introduction 


Plasma  propulsion  systems  are  more  advantageous  than  chemical  propulsion  systems  ac¬ 
cording  to  their  higher  specific  impulse  rates  and  are  more  advantageous  than  ion  engines 
based  on  their  less  complex  design  features.  Plasma  propulsion  systems  in  form  of  arcjets 
are  built  in  the  range  of  about  0.1  -  0.2  N  at  a  power  of  about  1  kW  as  station  keeping 
devices  and  in  a  thrust  range  of  about  1  -  2  TV  at  10  to  30  kW  for  drag  compensation 
(means)  and  orbit  maneuvering.  With  the  advent  of  onboard  nuclear  power  plants  of 
1  OOArW  electric  and  more  high  power  arcjets  and  even  MPD-thrusters  in  the  range  of  up 
into  the  MW  -range  could  offer  substantial  cost  reduction  or  payload  increase  in  space 
transportation  and  development. 

The  plasma  thrusters  which  are  considered  here  consist  of  coaxial  electrode  systems 
which  support  an  electrical  discharge  within  an  axial  gasious  propellant  stream.  The 
Joule  heated  gas  or  plasma  is  accelerated  by  thermal  expansion  and  by  self  induced  ele- 
tromagnetic  forces  and  depending  an  wether  one  or  the  other  accelerating  mechan.  :m 
dominates,  one  may  speak  from  arcjet-  or  MPD-  thrusters.  Indeed  the  most  advan¬ 
ced  plasma  thruster  and  almost  ready  for  space  flight  is  the  radiatively  cooled  1  kW 
[1,  2,  3,  4]  and  the  15  to  30  kW  arcjet  with  hydrazin  or  amonia  as  propellant  [5,  6,  7], 
the  less  advanced  thrusters  are  the  high  power  arcjets  while  MPD  thrusters  are  cur¬ 
rently  investigated  up  into  the  several  hundred  kilowatt  range  in  steady  state  mode  [7] 
and  within  the  megawatt  range  in  pulsed  mode  [8,  9].  Despite  of  their  simple  designs 
the  physical  proplems  involved  are  complex  and  by  no  means  fully  understood. 

The  objective  of  this  work  is  to  locate  and  to  investigate  the  still  unknown  effects  within 
the  flow  discharge  region  concerning  the  heating,  the  acceleration,  the  loss  mechanisms 
and  the  stability  behavior  of  the  arc  discharge  as  well  as  the  electode  attachment  pro¬ 
blematic  within  those  plasma  thrusters. 

Several  numerical  model  calculations  have  been  developed  which  allow  to  determine  the 
flow  arc  discharge  region  and  the  performance  (thrust specific  impulse)  of  a  nozzle  type 
and  a  cylindrical  MPD- thruster  as  afunktion  of  mass  flow  rate  and  electric  current.They 
are  described  in  chapter  2.  In  chapter  3  a  new  explanation  of  the  ”  Onset  "  -  Pheno¬ 
menon  (plasmastability  effect  in  MPD  thruster)  is  presented  which  is  based  on  a  "run 
away  Joule  heating  "  mechanism  of  the  current  carrying  plasma  channel.  The  results 
of  this  new  theory  are  compared  with  experimental  observations  and  are  in  very  good 
agreement. 

In  chapter  4  electrode  erosion  effects  were  investigated  experimentally  and  qualitatively 
explained  by  an  improved  spot  theory. 
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Figure  1.1:  Different  MPD-Thruster  Devices  a)  The  Nozzle  Type  Thruster  b)  The  Cy 
lindrical  Type  Thruster 


2.  MPD  -  Code  Development 


Based  on  the  two  different  self  magnetic  MPD  thruster  designs  a)  the  nozzle  type  and 
b)  the  zylindrical  type  thruster  as  shown  in  Fig.  1.1  a),  b);  two  different  MPD  -  codes 
have  been  developed  and  their  results  compared  with  experimental  ones. 

The  main  difficulties  for  developing  a  correct  physical  model  and  a  numerical  simulation 
are  the  following:  no  thermal  equilibrium  of  the  plasma,  real  gas  effects,  only  partly 
known  electrode  effects,  a  not  negligible  friction  due  to  the  low  Reynold’s  numbers  cau¬ 
sed  by  low  densities  and  high  temperatures,  a  magnetic  Reynold’s  number  too  high  to 
separate  flow  and  discharge,  etc.  These  are  only  the  most  important  physical  problems; 
added  to  them  are  the  numerical  difficulties;  different  sets  of  partial  differential  equa¬ 
tion,  namely  elliptical  and  hyperbolic  ones  that  have  to  be  coupled,  nonlinearities,  steep 
gradients  especially  at  the  electrodes,  difficult  geometries,  etc.  This  list  is  by  no  means 
a  complete  one. 

Due  to  this  complexity,  great  simplifications  have  to  be  made  for  all  codes  developed  to 
investigate  the  MPD  thrusters,  so  that  they  solve  at  least  parts  of  the  whole  problem. 

For  investigating  the  flow  and  discharge  fields  and  their  interaction  and  effects  on  electric 
power,  thrust  and  efficiency  of  the  MPD  thruster,  numerical  codes  have  been  developed. 
These  codes  are  one-dimensional  [10],  quasi-two-dimensional  [11,  12]  and  simplified  two- 
dimensional  ones  [1 3,  14,  15,  16]. 


2.1  The  Nozzle  Type  MPD  Thruster 


In  order  to  optimize  the  constrictor  size  to  a  given  power  and  mass  flow  rate,  three  nozzle 
type  thruster  configurations  (Fig  1.1  a)  have  been  investigated  experimentally  in  a  wide 
power  range  and  calculated  by  means  of  a  semi  two  dimensional  MPD  -  code.  This 
code  was  developed  to  predict  the  overall  performance  of  these  thrusters  and  to  prove 
the  validity  of  this  model  approach  by  comparing  the  computed  with  the  experimental 
results. 

In  this  approach,  the  quasineutral  part  of  the  discharge  is  modeled,  electrode  effects  are 
not  included.  For  this  purpose,  independent  codes  for  the  discharge,  the  flow  field  and 
the  electron  temperature  distribution  have  been  installed  and  then  solved  iteratively  (18). 
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The  relatively  low  magnetic  Reynold’s  number  in  the  region  of  magnetic  acceleration 
justifies  these  decoupled  calculations. 

The  basic  assumptions  for  all  codes  are  steady  state  conditions,  rotational  symmetry 
and  a  fully  singly  ionized,  quasi-neutral  plasma. 


2.1.1  Flow  field  code 

The  flow  field  is  taken  as  a  one-dimensional,  frictionless  expansion  flow,  assuming  adiaba¬ 
tic  behavior  for  the  ions  and  isothermal  for  the  electrones  within  the  entire  Nozzle.  One 
assumes  a  fully  ionized  plasma  and  distinguishes  between  electron  and  ion  temperatures 
and  assumes  a  frozen  flow  within  the  thruster.  Outside  the  thruster  both  components 
are  assumed  to  he  adiabatic.  The  results  of  the  codes  are  compared  with  data  from 
experiments  performed  within  a  vacuum  tank  which  always  has  a  residual  gas  pressure. 
Therefore,  a  free  stream  boundary  is  implemented  corresponding  to  experiments  with  a 
similar  thruster  [19].  Outside  the  free  stream  boundary,  see  Fig.  2.1,  the  ambient  gas 
mass  density  is,  according  to  the  experiments,  set  to  10~6  corresponding  to  a  pressure 
of  0.1  Pa.  At  the  free  stream  boundary  a  linerar  transition  of  the  mass  density  over  a  di¬ 
stance  of  30mm  is  implemented  to  avoid  numerical  problems  by  the  solution  of  equation 
2.12,  caused  by  discontinuities. 

The  assumptions  for  the  flow  field  code  are:  a  one-dimensional  frozen  flow,  an  ambipolar 
expansion  flow  and  no  friction. 

The  basic  equations  are: 


•  Continuity  Equation: 

•  Bernoulli  Equation: 

•  Equation  of  State: 

-  Within  the  thruster: 

-  Outside  the  thruster: 


gvA  =  m 


v  dv  +  -dp  =  0 
0 


(2.1) 

(2.2) 


(2.3) 

(2.4) 


pc  and  qc  are  the  chamber  pressure  and  chamber  density,  respectively;  they  can  be 
calculated  from  the  nozzle  throat  condition. 
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Within  the  nozzle  throat  one  assumes  that  all  the  ohmic  heat  input  per  unit  arc  length 
is  dissipated  radially  by  heat  conduction.  This  assumption  is  quite  reasonable,  since 
the  axial  temperature  gradient  is  much  smaller  than  the  radial  one  and  since  any  heat 
input  at  the  smallest  cross  sectional  area  cannot  increase  the  velocity  beyond  the  speed 
of  sound.  The  radiation  loss  can  be  neglected,  because  within  an  MPD  thruster  the 
chamber  and  throat  pressures  are  low. 

A  detailed  evaluation  of  this  energy  balance  within  the  nozzle  throat,  made  in  [18],  yields 
the  following  equations  for  calculating  the  electron  temperature  within  the  nozzle 

r.  =  /(-)  (2.5) 

where  /  is  a  slowly  varying  function  depending  on  the  current  distribution. 

The  equation  of  state  together  with  the  continuity  and  Bernoulli  equations  yield  the 
following  expansion  equations: 


-  inside  the  thruster: 


^(^(?)  ♦■-(in -&)■-•  “ 

-  outside  the  thruster: 

*+AE(‘ -(£)”') -(a)'-  ,2I) 

where  veo,  peo  and  are  the  values  at  the  thruster  nozzle  exit. 

The  results  of  this  flow  field  code  arc  the  axial  velocity,  pressure  and  density  distribution. 

After  the  first  step  of  iteration  of  all  codes,  when  the  current  density  field  is  known,  the 
pressure  rise  due  to  the  pinch  effect  has  been  taken  into  account.  This  results  in  the 
pressure  profile: 

P(r)  =  p(R)  —  (' (J  x  B)r  dr  (2.8) 

JR 

Together  with  the  continuity  equation, one  obtains  the  density  profile. 


2.1.2  Discharge  code 

In  order  to  calculate  the  current  distribution  of  the  arc  discharge,  a  two-dimensional 
computer  code  has  been  developed. 

The  additional  assumptions  for  the  discharge  code  is  the  neglection  of  azimuthal  current. 
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(2.9) 


The  basic  equations  for  the  discharge  is  the  extended  Ohm’s  law  for  plasmas: 

J=<j(E  +  vx  B)  -  ^(/*  B) 

Rewriting  the  Ohm’s  law  by  means  of  Maxwell’s  equations,  one  obtains  a  vector  equation 
for  the  magnetic  induction  vector  B  in  the  form 

0  =  —(V  x  (-V  x  B)  -  (V  x  (w  x  B)) 

Vo  <* 


+  — (V  x  (0(V  x  B)  X  §)) 

Vo 

(2.10) 

0=~-~ 
dbg  ene 

(2.11) 

From  this  equation  follows  with  a  stream  function  =  r  Be  and  with  respect  to  the 
rotational  symmetry  and  the  zero  azimuthal  current  the  elliptical,  partial  differential 
equation  of  2nd  order 

d21>  dH 
dr 2  dz2 

dV  ( 1  1  da  dd 

-~rr  -  +  -T5 - a-  +  <?VoVr 

dr  \r  a  dr  r  dz 

d'P  (Ida  a  dfl  2<7/?’F 

dz  \<7<9c  +  r  dr  r 2  +  ^ 

.  ( dvT  dv,  v,\ 

(2-,2> 

The  function  >F( r,z )  =  const  represents  now  a  current  contour  line,  since  B  =  Be  is 
proportional  to  where  7(r)  is  the  electric  current  carried  through  a  cross  sectional 
area  of  nr2.  The  proper  boundary  conditions  for  <F  follow  from  the  geometry  of  the 
thruster  walls  and  electrodes.  At  the  insulator  inside  the  thruster  and  at  the  inflow 
boundary  <F  is  set  to  —  For  the  electrodes  one  assumes  that  the  electric  field  is 
normal  to  the  surfaces,  therefore;  E  ■  t  —  0,  where  t  is  the  tangential  vector  of  the 
electrodes  surfaces  in  r-z  plain.  At  the  other  boundary  sections  'F  is  set  to  0. 


2.1.3  Electron  temperature  code 

The  electron  temperature  has  a  strong  effect  on  the  electrical  conductivity  and  the  elec¬ 
tron  density,  which  again  influences  the  discharge  pattern.  Therefore,  a  two-dimensional 
code  for  the  electron  temperature  distribution,  corresponding  to  the  two-dimensional 
discharge  code,  has  been  written. 
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The  electron  temperature  distribution  is  determined  by  the  energy  equation  for  the 
electron  component. 

-  =  — AeV  •  Te  -  §-/Vr£  (2.13) 

<7  2  e 

The  energy  input  due  to  ohmic  heating  equals  the  sum  of  losses  given  by  the  heat  flux 
due  to  conduction  and  convection.  The  thermal  conductivity  Ae  depend  on  the  elektron 
temperature: 


Ae  =  2.08 


•JmlQc 

where  Qce  is  the  Gvosdover  cross  section  [19] 

.  2 


\"T$ 


=  *-(- 


4  \47TCo kTe 


+  1447T' 


4T?k 


nce6(zf(zi  +  1))> 


(2.14) 


(2.15) 


With  respect  to  the  rotational  symmetry,  equations  2.13  -  2.15  result  in  the  following 
elliptical,  partial  differential  equation  of  2nd  order: 


(2.16) 


The  boundary  condition  for  Tt  along  the  free  stream  boundary  is  set  to  a  value  of  10000 
K  in  accordance  with  the  measurements  [19].  Outside  the  free  stream  boundary  the 
electron  temperature  is  set  to  a  constant  value  of  7000  K  [19].  To  avoid  a  discontinuity 
at  the  free  stream  boundary,  a  linear  transition  in  radial  direction  is  used.  By  using  a 
value  of  10000  K  for  the  anode  and  the  insulator  one  obtains  reasonable  results  for  the 
electron  temperature  inside  the  nozzle.  The  cathode  is  treated  as  a  thermal  insulator; 
therefore  VTe  ■  n  =  0,  where  n  is  the  normal  vector  of  the  cathode  surface.  Due  to  the 
axial  symmetry,  ^  =  0  on  the  axis. 


2.1.4  Solution  of  the  Equations 

The  three  models  described  above  were  solved  in  the  following  manner: 
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For  a  given  flow  field,  the  discharge  equation  and  hence  the  electron  temperature  are 
determined.  With  these  results,  in  the  next  iteration  step,  the  flow  field  equations  are 
calculated.  This  new  flow  field  and  the  new  distribution  of  the  electron  temperature 
have  been  taken  to  calculate  the  discharge  equation,  and  so  on.  This  iteration  steps 
are  repeated  until  the  calculated  datas,  such  as  electron  temperature  distribution  and 
current  distribution  have  reached  numerical  equilibrium.  Usually  eight  iterations  are 
necessary  to  get  good  results. 

The  expansion  relations  2.6  and  2.7  have  been  solved  with  the  Newton  method.  The 
two  nonlinear,  elliptical,  partial  differential  equations  2.12  and  2.17  are  solved  with  the 
finite  differential  method.  Figure  2.1  shows  the  thruster  with  the  free  stream  boundary 
which  devides  the  field  into  two  areas.  Area  1  and  2  is  used  as  calculation  area  for  "P, 
whereas  only  area  1  is  used  for  Tt.  These  areas  are  discretizised  with  a  aquidistant  finite 
differential  grid  of  1mm  height  (radial)  and  2mm  width  (axial).  The  partial  derivatives 
of  1st  and  2nd  order  are  replaced  by  the  following  differential  quotients  of  any  selective 
function  f(x) 


df  _  /x+i  ~  fx—i 

dx  2Ax 

d2l  =  /»+! -/.  +  /*-! 

dx 2  Ax 2 


(2.17) 

(2.18) 


With  this  discretization,  one  obtains  a  nonlinear  system  of  equations,  which  is  iteratively 
solved  with  a  modified  Gauss-Seidel  algorithm. 


2.1.5  Discussion  of  Results 

A  main  goal  for  developing  this  numerical  code  is  to  obtain  better  insight  into  the 
current  distribution  within  a  continuously  running  nozzle  type  MPD  thruster.  With 
that  developed  code,  the  three  geometries  as  shown  in  Fig.  1.1a  have  been  calculated. 

With  quasi-steady  pulsed  thrusters  it  is  possible  to  measure  the  current  contour  lines, 
but  it  is  not  yet  possible  to  compare  these  with  experimental  data  on  the  continuously 
running  thruster  because  of  the  high  heat  load  of  the  probes.  This  means  there  is  a  need 
for  a  calculation  code  for  the  steady  state  thruster  development. 

The  magnetic  Reynold’s  number  indicates  if  the  discharge  influences  the  flow  field  si¬ 
gnificantly.  Therefore  to  evaluate  the  only  weak  coupling  of  flow  and  discharge  within 
our  iterative  solution  of  the  three  codes,  the  magnetic  Reynold’s  number  Rm  =  fiocrvL 
has  been  calculated.  In  Fig.  2.2  its  distribution  for  the  three  throat  geometries  within 
the  discharge  region  has  been  plotted,  where  as  reference  length  the  throat  diameter  was 
chosen.  Since  the  main  influence  of  the  electromagnetic  forces  on  the  flow  field  occurs  in 
the  supersonic  part,  where  the  magnetic  Reynold’s  number  is  fairly  low,  the  decoupling 
of  the  flow  field  and  the  electromagnetic  field  model  seems  justified. 
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A  comparison  between  calculated  and  measured  current  contour  lines  has  been  shown 
earlier  [12].  The  comparison  shows  fairly  good  agreement  with  the  DT2  thruster  con¬ 
figuration.  In  Fig.  2.3  the  calculated  current  contour  lines  are  shown  for  the  different 
throat  geometries  at  2000  A  with  a  mass  flow  of  0.8  g/s  argon. 

The  electron  temperature  distribution  within  and  outside  the  thruster  at  a  current  level 
of  2  kA  is  shown  in  Fig.  2.4.  The  maximums  of  the  electron  temperatur  in  the  symmetry 
axis  are  in  the  middle  of  the  throat,  because  of  the  high  current  density  there.  The 
maximum  electron  temperatures  are:  34  600  K  for  DT5,  37  700  K  for  DT2  and  51 
100  K  for  DT3.  The  axial  temperature  profile  is  compared  with  experiments  [19]  in 
Fig.  2.5.  The  coincidence  is  fairly  good.  The  reason  for  the  lower  calculated  temperatures 
upstream  is  that  the  temperature  boundary  condition  on  the  upstream  flow  boundary 
of  10000  K  at  a  distance  of  185  mm  from  the  thruster  end  is  too  low  compared  with  the 
experiment. 

The  potential  lines  within  the  discharge  at  2  kA  are  plotted  in  Fig.  2.6.  An  integration 
across  these  lines  yields  the  discharge  voltage. 

This  calculated  discharge  voltage  is  compared  with  the  measured  voltage  in  the  current 
range  between  1  and  5  kA  in  Fig.  2.7. 

With  the  measurements,  in  contrast  to  the  calculation,  the  electrode  fall  voltages  are 
included.  In  the  case  of  the  cathode  this  voltage  drop  equals  about  the  work  function 
of  the  cathode  material,  which  means  in  the  case  of  thoriated  tungsten  between  1  and 
3  volts,  and  is  nearly  independent  of  the  discharge  current.  The  anode  fall  voltage  de¬ 
pends  strongly  on  the  current;  starting  with  several  negative  volts  at  currents  up  to  about 
3  kA,  the  anode  voltage  rises  steeply  in  the  onset  region  [10].  This  explains  the  different 
behavior  of  the  two  curves  at  high  currents  and  the  fairly  good  agreement  for  currents 
between  1.5  and  3  kA.  For  DT3  it  was  not  possible  to  measure  with  currents  over  1,6  kA. 


According  to  the  flow  conditions  and  the  known  electromagnetic  force  configuration,  the 
thrust  can  be  calculated.  The  thrust  of  an  MPD  thruster  is  the  sum  of  all  gas  dynamic 
surface  forces  and  the  electromagnetic  volume  forces.  Hence  it  is 

f  =  fA  (evv  +  p)-AA-  Jvjx  B  dV  (2.19) 

where  A,  represents  the  surface  of  all  internal  walls  and  V  is  the  current  carrying  volume. 
The  thermal  loss  to  the  wall  of  the  thruster  were  neglected  in  this  equation.  Using  the 
well-known  method  to  match  the  gasdynamic  effect  by  using  a  thrust  coefficient  cr, 
where  ct  is  defined  as 

c’=^A.LX,'ss+f)iX  (2'M) 

the  absolute  value  of  the  thrust  yields  then 

I  f  |=  crPcA,+  |  Jv  j  x  §  dV  |  (2.21) 
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As  it  is  shown  in  Fig.  2.8,  one  obtains  for  all  three  devices  DT2,  DT3  and  DT5  a  very 
good  fit  between  measured  and  calculated  thrust  by  taking  a  thrust  coefficient  of  ct  =  1. 


To  improve  the  flow  data  at  the  upper  flow  boundary,  which  were  calculated  from  the 
nozzle  throat  model,  the  calculated  chamber  pressure  is  compared  with  measurements 
carried  out  with  DT2-IRS  (Fig.  2.9).  The  very  low  deviation  which  is  always  less  than 
10  per  cent  is  an  indication  of  the  validity  of  this  model. 

Since  the  pinch  effect  is  a  two-dimensional  phenomenon  and  is  applied  here  numerically 
on  a  one-dimensional  flow  model,  numerical  instabilities  occur  with  currents  higher  than 
3.5  kA  at  a  mass  flow  of  0.8  g/s.  These  difficulties  can  only  be  overcome  with  a  two- 
dimensional  flow  model.  The  results  for  4  kA  and  higher  were  therefore  calculated 
without  accounting  for  the  pinch  effect.  To  evaluate  the  influence  of  the  pinch  effect 
on  the  discharge,  a  calculation  was  carried  out  with  and  without  respect  to  the  pinch 
effect  at  a  current  of  3  kA  and  a  mass  flow  of  0.8  g/s  argon  for  comparison  with  the 
results  of  the  codes  described  above.  No  significant  influence  could  be  detected  on  the 
electron  temperature  distribution,  the  discharge  voltage  or  the  thrust.  The  influence  on 
the  current  distribution  is  demonstrated  in  Fig.  2.10.  The  main  difference  is  that  with 
respect  to  the  pinch  effect  the  current  is  at  the  anode  more  concentrated  at  the  nozzle 
end. 


2.1.6  Conclusions 

The  presented  numerical  analysis  has  been  applied  to  a  rather  complex  nozzle  type 
geometry.  Despite  the  shortcomings  of  the  code,  namely  the  decoupled  computation 
of  the  two  dimensional  electromagnetic  and  the  one  dimensional  flow  equations  yields 
valuable  results: 


•  Calculated  integral  values,  such  as  thrust,  discharge  voltage  and  discharge  chamber 
pressure  fit  the  measured  ones  very  well. 

•  The  current  patterns  of  the  discharge  and  the  electron  temperature  coincide  with 
the  measured  data  fairly  well. 


Further  steps  to  improve  the  method  will  be  a  refinement  of  the  flow  field  code  to  two- 
dimensionality  and  the  inclusion  of  chemical  reactions  to  overcome  the  restriction  of  a 
fully  singly  ionized  gas  assumption. 
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2.2  The  Cylindrical  MPD  Thruster 


A  new  fully  two  -  dimensional  axissymmetric  code  has  been  developed  in  which  the 
gas  dynamic  equations  are  solved  by  the  finite  volume  code  EULFLEX  [21]  The  elec¬ 
tromagnetic  discharge  equation  and  the  electron  energy  equation  are  transformed  into 
curvilinear  coordinates  and  solved  under  steady  state  conditions  by  a  modified  Gauss- 
Seidel  algorithm. 

In  order  to  weigh  the  suitability  of  that  code,  it  is  used  to  investigate  the  cylindrical 
MPD  thruster  (Fig.  1.1  b),  which  is  investigated  experimentally  at  the  IRS  in  a  steady 
state  as  well  as  in  a  quasistationary  pulsed  mode. 

All  geometrical  configuration  parameters  on  which  the  computation  is  based  on  were 
chosen  in  accordance  with  the  experimental  device.  A  transformation  to  curvilinaear 
coordinates  allows  a  calculation  of  a  complex,  two-dimensional  geometry.  Fig.  2.11 
shows  a  grid  of  those  calculations  for  the  cylindrical  thruster  .  The  propellant  used  in 
the  calculation  code  is  argon. 

The  two-dimensional  numerical  code  is  based  on  a  finite  volume  method  for  flow  calculations. [21] 
and  the  Gauss- Seidel  method  for  the  electromagnetic  discharge^  1,  21]  and  the  electron 
energy  equation[23,  24]. 

In  addition  to  the  geometry  of  the  thruster,  the  code  requires  the  mass  flow  rate  and 
the  electric  current  as  input  data,  at  the  starting  point  the  heavy  particles  temperature 
and  Mach  number,  and  as  an  input  from  the  experimental  investigations  the  boundery 
values  of  the  electron  temperature.  This  yields  the  following  results:  flow,  electron-  and 
heavy  particles  temperature  distribution,  pressure  distribution,  velocity,  Mach  number 
distribution,  magnetic  field,  current  density  distribution  and  the  electromagnetic  and 
electrothermal  thrust. 


In  order  to  calculate  the  current  distribution  of  the  arc  discharge,  a  two-dimensional 
computer  code  has  been  developed.  Additional  assumptions  for  the  discharge  code  are 
the  neglection  of  azimuthal  current. 


The  basic  equation  for  the  discharge  is  the  extended  Ohm’s  law  for  plasmas: 

J  =  a(E  +  v  x  B)  -  ^(/  x  B)  (2.22) 

Rewriting  the  Ohm’s  law  by  means  of  Maxwell’s  equations,  one  obtains  a  vector  equation 
for  the  magnetic  induction  vector  B  in  the  form 


with 


0  =  —(V  x  (iv  x  B)  -  (V  x  (vx  B)) 
Po  <r 

+— (V  x  (/?(V  x  B)  x  B )) 
f*  o 


1 

ene 


(2.23) 

(2.24) 
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With  the  Jacobian  J  of  the  transformation  from  cylindrical  to  the  curvilinear  coordinates: 


J  “  d!dr_  _  dzdr_  (2.25) 

di  dv  dyj  di 

where  the  matrices  are  formed  in  terms  of  the  derivatives  of  the  cylindrical  coordinates 
as  follows: 


dr 

.dr 

6 

-  J-k-  , 
dr) 

^  =  ~Jdi 

,dz 

wdz 

6 

=  -J~-  , 

dij 

r,T  =  Jdl 

where  the  convention  =  |£,  etc.  is  used. 


(2.26) 


The  equation  (2.23)  follows  with  a  stream  function  =  r  Bg  the  elliptical,  partial 
differential  equation  of  2nd  order 


7&T  (£  +  62)  +  +  6»?r)  +  TO-  +  *>') 


d2* 


0(0 

+ 


'dtdjf 


Or? 


d* 

d( 


2(70* 
R 2 


-  (7  HO  U 


(2.27) 
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with  U  and  V  as  the  contravariant  velocities 


U  -  £;U+  (rv 

V  =  i}zu  +  i)Tv  (2.28) 


The  function  <I»(r,  z)  =  const  now  represents  a  current  contour  line,  since  B  —  Be  is 
proportional  to  where  /(r)  is  the  electric  current  carried  through  a  cross  sectional 
area  of  7rr2.  The  proper  boundary  conditions  for  follow  from  the  geometry  of  the 
thruster  walls  and  electrodes.  At  the  insulator  inside  the  thruster  and  at  the  inflow 
boundary  'f  is  set  to  —  fj/.  For  the  electrodes  the  electric  field  is  assumed  to  be  normal 
to  the  surfaces  E  ■  f  =  0.  In  accordance  with  the  grid  in  Fig.  2.1 1  it  is  equal  to  =  0, 
and  it  follows  to: 


di) 


TT  (SrVr  -  +  -  (f zVr  ~  i^z) 

n  o 


di  R 

At  the  other  boundary  sections  and  at  the  symmetry  axis  'P  is  set  to  0. 


(2.29) 


The  electron  temperature  has  a  strong  effect  on  the  electrical  and  thermal  conductivity 
and  on  the  electron  density,  which  again  influences  the  discharge  pattern.  Therefore, 
a  two-dimensional  code  for  the  electron  temperature  distribution,  corresponding  to  the 
two-dimensional  discharge  code,  was  written. 


The  electron  temperature  distribution  is  determined  by  the  energy  equation  for  the 
electron  component. 


1  k 

V(A  «vr.)  +  = 


(T)  -  Trot) 

V 


Mon 


(2.30) 


The  first  term  on  the  left  hand  side  represents  the  heat  flux  of  the  electron  gas  due  to 
condution,  and  the  second  term  gives  the  convective  change  of  the  electron  gas  energy. 
The  energy  input  due  to  ohmic  heating  is  represented  by  the  first  term  on  the  right  hand 
side.  The  sum  of  losses  due  to  the  energy  transfer  from  the  electron  gas  to  the  heavy 
particles  gas  is  calculated  by  the  second  term  on  the  right  hand  side.  The  reaction  losses 
are  given  by  the  last  term,  where  e,on  is  given  by  Unsold. (25) 


The  heat  transfer  coefficient  aeu  and  the  thermal  conductivity  A„  depend  on  the  electron 
temperature  [23,  28],  where  the  heat  transfer  coefficient  follows  by 


8>/2 

\T* 


kT't 


771  e  -f  771  y 


and  the  thermal  conductivity  by 


15  [W  nek7Te _ 

8  V  2  E^=e)  nvQtvy/mtvkTe v 


(2-31) 


(2.32) 
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where  for  Qec  and  for  Qeu  with  u  =  ions  the  Gvosdover  cross  section[26]  is  taken,  and 
with  i/=atoms  the  Ramsauer  cross  section  is  applied. [27] 


With  respect  to  the  rotational  symmetry,  equation  2.30  results  in  the  following  elliptical, 
partial  differential  equation  of  2nd  order: 


*  d? 


+  20 


d2Te 


+  7' 


a2T. 


d£dri  di] 


+ 


J2 

4Aecr//2 


^nen„a„(7'I/  -  Tc)  + 


3  nek 

~2f 


j  (  3ne  ret  .  ret  \ 

- c.on  •  M—gT—  +  v— 5 -  - 

°  \  d(,  dri  ) 

nekTe 

—j—  (Wt  -  ffU,  +  2(V n  -  r„t;e)  - 


( 


dTe  8TC  dTt 
Ur^-l,r^  +  VZ^'VZ" 


OTA 

dZ). 


(2.33) 


Here  (,  and  q  indicate  the  partial  differentiation  to  the  curvilinear  £-  and  ^-directions. 
a,  0,  7,  er,  r  are  metric  terms  and  J  is  the  Jacobian  of  the  transformation. 


At  the  outflow  boundary  Te  is  set  to  a  value  in  accordance  with  measurements[19].  Inside, 
the  solid  bodies  of  the  thruster  are  treated  as  thermal  insulators;  therefore  VTC  •  n  =  0. 
where  n  is  the  normal  vector  of  the  surfaces.  Due  to  the  axial  symmetry,  =  0  on  the 
axis.  At  the  inflow  boundary  T,  is  set  to  a  constant  value  of  7000  A'  in  accordance  with 
the  measurements. 


For  the  description  of  the  two-dimensional,  axisymmetric  flow,  the  following  nonlinear 
hyperbolic  system  of  differential  equations  with  curvilinear,  cylindrical  coordinates  is 
used. 

§1  Q  +  ^  F(q)  +  ^  G(q)  +  H(<)  =  0  (2.34) 

The  indices  t,  <,  r,  indicate  the  partial  differentiation  with  respect  to  time  and  to  the  £- 
and  ^-directions. 


The  first  three  terms  are  used  in  the  usual  fluid  dynamic  manner,  where  q  is  the  flow 
variables  vector  and  F  and  G  are  the  spatial  derivatives  vectors. [21]  The  source  vector 
//  is  developed  to; 


’ 

pv/R 

* 

puv/R+^L 

)  +  j>7 

ndV  ,rdU  ,8V\ 
L  9r,  T  i  di,  ^  d(  > 

/  + 

L  (p  +  e)v/R  +  Vq,  -  £„  n„nea^(T -Tt)  J 

(2.35) 


The  first  terms  in  the  source  vector  transform  the  plane  two-dimensional-  into  a  cylindri¬ 
cal  calculation.  The  first  additional  terms  in  the  impulse  equations  represent  the  stress 
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tensor  of  the  plasma  flow,  where  the  viscosity  coefficient  //  is  given  by: 

__  3  fit  _ mtntkTt _  fo^A) 

11  ~  8  V  2  Y  E,(=0  nuQi„\/mlukTlu 

The  last  terms  in  the  impulse  equations  represent  the  j  x  B  forces  from  the  electroma¬ 
gnetic  field.  In  the  energy  equation  the  second  source  vector  term  qs  represents  the  heat 
flux  vector  due  to  the  heavy  particles  of  the  MPD  flow,  which  is  given  by  [28]: 

£  =  L-A4-V»<  +  V7i)  (2.37) 

l  £  Til 

Here  the  heat  conductivity  coefficients  are  given  by: 

t  f*  nik*Tl  10 

"  8V2 

The  transfer  of  Joules  heat  is  contained  in  the  temperature  exchange  between  the  electron 
and  the  heavy  particle  temperature,  which  is  represented  by  Y.nuncacu(Tt  —  T)  in  the 
equation  of  energy  conservation. 

The  transport  coefficients  are  derived  more  detail  in  [28].  The  cross  sections  Qi„  for  all 
coefficients  arc  taken  from  the  references. [26,  27]. 


2.2.1  Solution  of  the  Equations 

The  three  models  described  in  the  subsection  above  are  solved  in  the  following  manner: 

For  a  given  flow  field,  the  discharge  equation  and  hence  the  electron  temperature  is 
determined.  With  these  results  the  flow  field  equations  are  calculated.  This  new  flow 
field  and  the  new  distribution  of  the  electron  temperature  are  used  to  calculate  the 
discharge  equation,  and  so  on  (see  Fig.  2.12).  These  iteration  steps  are  icpeated  until 
the  calculated  datas,  such  as  electron  temperature  distribution,  current  distribution  and 
the  heavy  particles  flow  field,  have  reached  numerical  equilibrium. 

The  extended  Euler  eqations  which  determine  the  heavy  particles  flow,  are  solved  by 
the  finite  volume  method  EUFLEX.[21]  The  two  nonlinear,  elliptical,  partial  differential 
equations  2.27  and  2.33  are  solved  with  a  Gauss-Seidel  finite  difference  method. 


2.2.2  Numerical  Results 

The  magnetic  Reynold’s  number  in  Fig.  2.13  is  based  on  the  distance  between  the 
electrodes.  It  indicates  if  the  induced  magnetic  field  influences  the  flow  field  significantly. 
In  Fig.  2.13  one  can  see  that  the  magnetic  Reynold’s  number  increases  only  at  the  end 
of  the  cathode  to  1.1;  in  the  greater  part  of  the  discharge  area  it  is  considerably  lower 
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then  1.  Therefore  the  weakly  coupled  calculation  of  flow  and  discharge  as  performed 
here  leads  to  fairly  good  results. 

The  flow  inlet  boundary  conditions  were  iterated  to  coincide  with  an  experimentally 
obtained  cold  gas  thrust  of  1  N  at  a  mass  flow  rate  of  2 g/s.  For  a  given  current  of  6kA 
the  computation  yields  the  current  contour  lines,  as  illustrated  in  Fig.  2.14. 

This  calculated  current  distribution  corresponds  to  that  of  a  continuosly  running  thruster 
in  which  the  cathode  is  hot  glowing  and  emits  sufficient  electrons  along  its  entire  length. 

In  case  of  pulsed  thruster  with  cold  cathode  the  arc  preferably  touches  the  cathode 
root,  therefore  the  current  contour  distribution  is  not  comparable  with  the  case  of  the 
continuing  mode.  Unfortunately  it  is  only  possible  with  quasi-steady  pulsed  thrusters 
to  measure  the  current  contour  lines.  Due  to  the  high  heat  load  to  the  probe  it  is  not 
yet  possible  to  compare  the  calculated  with  the  experimental  data  on  the  continuously 
running  thruster.  In  order  to  verify  the  calculation  procedure,  an  experiment  with  a 
withdrawn  cathode  was  carried  out,  since  here  the  arc  attachment  is  closer  to  the  top  of 
the  cathode  also  in  case  of  a  quasi  steady  state  mode.  In  that  case  it  means  there  is  a 
need  for  a  calculation  and  measurement  with  a  short  cathode,  because  in  that  position 
the  current  has  hardly  any  space  to  move,  and  the  cathode  attachement  of  the  arc 
during  a  measurement  is  well  defined.  This  promises  a  better  correspondence  with  the 
calculations. 

Due  to  a  limitation  of  the  pulsed  power  supply  a  comparison  betwer  i.  ensured  andcalcu- 
lated  current  contour  lines  has  been  done  only  at  3000  /I  for  a  c\,^hodc  position  of  30  mm 
(see  Fig  2.15).  For  this  condition,  however,  the  <omparision  is  not  quite  appropriate, 
since  the  difference  between  pulsed  and  steady  state  thrusters  is  much  more  distinct 
for  the  cylinder  type  thruster  than  for  the  nozzle  type  thruster.  In  the  experiment  the 
current  touches  mainly  the  upstream  area  of  the  cathode  whereas  the.  calculation  results 
in  a  stronger  current  concentration  at  the  tip  of  the  cathode.  These  different  current 
distributions  on  the  cathode  cause  a  more  downstream  current  on  the  anode  for  the 
calculated  case  than  in  the  experiment. 
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Figure  2.2:  Magnetic  Reynold’s  Number  Distribution  of  the  Three  Different  Throat 
Geometries  for  2  kA  and  a  Mass  Flow  of  0.8  g/s  Argon.  As  Characteristic  Length  the 
Nozzle  Diameter  was  Chosen. 
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Figure  2.3:  Calculated  Current  Contour  Lines  of  the  Three  Different  Throat  Geometries 
for  2  kA  and  a  Mass  Flow  of  0.8  g/s  Argon. 
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Figure  2.4:  Electron  Temperature  Distribution  of  the  Three  Different  Throat  Geometries 
for  2  kA  and  a  Mass  Flow  of  0.8  g/s. 
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Figure  2.6:  Potential  Lines  of  the  Three  Different  Throat  Geometries  for  2  kA  and  a 
Mass  Flow  of  0.8  g/s. 
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Figure  2.7:  Calculated  and  Measured  Discharge  Voltage 
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Figure  2.8:  Comparison  of  Calculated  and  Measured  Thrust  for  a  Mass  Flow  of  0,8  g/s 
Argon  cj  =  1 
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Figure  2.9:  Calculated  and  Measured  Arc  Chamber  Pressure 


Figure  2.10:  Density  Maps  a)  with  and  b)  without  the  Pinch  Effect  of  the  DT2  with 
Argon  Discharge  at  2  kA  and  a  mass  Flow  of  0.8g/s 


C  mruD 


Figure  2.13:  Magnetic  Reynold’s  Number  Distribution. 


Figure  2.14:  Calculated  Current  Contour  Lines. 
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Figure  2.15:  Measured  a)  and  Calculated  b)  Current  Contour  Lines. 


The  heavy  particle  (a)  and  the  electron  temperature  distribution  (b)  within  and  outside 
the  thruster  is  shown  in  Fig.  2.16.  The  maximum  temperature  value  of  the  heavy 
particles  and  of  the  electrons  is  near  the  tip  of  the  cathode,  but  there  is  also  a  remarkable 
increase  of  the  electron  temperature  at  the  upstream  side  of  the  anode. 

The  density  map  (Fig.  2.17)  demonstrates  the  expansion  flow,  where  the  relatively  low 
pinch  effect  has  been  taken  into  account. 

The  relation  of  the  velocity  distribution  to  the  temperatures  and  density  distribution  is 
shown  in  Fig.  2.18  as  a  vector  graph.  It  shows  a  high  increase  from  the  inflow  boundary 
downstream.  The  radial  velocity  components  inside  the  channel  are  quite  small,  which 
is  also  an  indication  for  a  relatively  low  pinch  effect  there. 

According  to  the  flow  conditions  and  the  known  electromagnetic  force  configuration,  the 
thrust  can  be  calculated.  The  thrust  of  an  MPD  thruster  is  the  sum  of  all  gas  dynamic 
surface  forces  and  the  electromagnetic  volume  forces.  Hence  it  is 

f  =  jf  (e  VV  +  p)  ■  dA  -  f  j  x  B  dVc  (2.39) 

where  ,4S  represents  the  surface  of  all  internal  walls  and  Vc  is  the  current  carrying  volume. 
For  the  case  presented  here  the  total  thrust  was  calculated  without  any  losses  as  14.4  N. 
This  includes  the  portion  of  5.1  N  electromagnetic  thrust  and  the  portion  of  1  N  cold 
gas  thrust. 


2.2.3  Comparison  with  experimental  results 

As  can  be  seen  in  Fig.  1.1b  the  anode  of  the  cylindrical  thruster  is  divided  into  three 
segments  of  30  mm  width.  The  current  in  each  of  the  segments  is  measured  seperately 
during  the  experiment. 


|  Percentual  Current  Fraction  into  Anode  Segments 

Discarge  Voltage 

Anode  1 

Anode  2 

Anode  3 

Theory 

49  A 

29  A 

22  A 

Experiment 

62  A 

25  A 

13  A 

mm 

Table  2.1:  Comparison  of  Numerical  and  Experimental  Results  for  2  g/s  Argon  at  6000 
A 


The  measured  data  of  table  2.1  agree  well  with  the  numerical  results  of  Fig  2.14.  For 
example,  the  experimentally  determined  current  percentage  are  62%  of  the  total  for  the 
downstream  anode  segment,  25%  for  the  middle,  and  13%  for  the  upstrem  section.  This 
must  be  compared  with  the  theoretical  percentage  of  49%,  29%  and  22%.  The  measured 
voltage  drop  was  18  V,  whereas  the  calculation  gave  a  somewhat  lower  figure  of  13  V. 
Since  the  numerical  model  does  not  include  electrode  fall  voltages,  the  resulting  number 
matches  the  experiment  quite  well. 
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The  code  described  in  this  paper  was  developed  and  tested  in  order  to  gain  a  more 
profound  understanding  of  the  fundamental  processes  occuring  in  MPD  thrusters  and 
in  order  to  end  up  with  reliable  design  criteria  and  predictions  of  thruster  performances 
under  various  conditions.  With  the  curvilinear  grid,  the  code  permits  a  more  appro¬ 
ximated  simulation  of  geometry  dependencies.  Thus,  better  datas  are  obtained  for  the 
calculation  results  with  the  cylindrical  thruster  design.  In  a  later  research  period  also 
the  nozzle  type  thruster  designs  will  be  investigated  with  the  curvilinear  code. 
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Figure  2.16:  a)  Heavy  Particles  and  b)  Electron  Temperature  Contours 
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Figure  2.18:  Velocity  Vector  Distribution. 
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3.  Plasmastability 


The  plasmastability  within  a  self  magnetic  MPD  -  thruster  has  been  or  is  being  investiga¬ 
ted  experimentally  and  theoretically  many  times.  One  of  the  most  crucial  limiting  effets 
of  these  thrusters  turns  out  of  to  be  the  so  called  "Onset”  phenomenon  [31,  32,  33,  34], 
a  plasmainstability  which  causes  a  severe  degradation  of  the  thruster  performance.  Se¬ 
veral  different  explanations  [35,  36,  37,  38,  39]  have  been  put  foreward  so  far  and  their 
analytical  results  agree  for  almost  all  quite  well  with  experimental  observations.  Here, 
another  explanation  will  be  presented  which  also  leads  to  exellent  agreement  with  the 
experimentally  found  onset  and  which  may  have  some  merits.  It  is  based  on  the  ”were- 
dependency”  of  the  electron  heatflux  (see  Appendix  A.l)  and  the  fact  that  with  incre¬ 
asing  uictc  an  overheating  of  the  discharge  channel  may  occure.  Such  a  process  which 
consequently  would  lead  to  an  unstable  "run  away  heating”  could  be  avoided  if  the 
channel  assumes  a  somewhat  curved  e.  g.  a  helical  configuration  thereby  leading  to  an 
increased  convection  cooling  effect  [magnetic  pumping  effect  [38]].  In  that  case  increased 
heating  inside  the  channel  would  be  balanced  by  an  increased  loss  mechanism. 

In  the  following  calculation  the  condition  at  which  such  an  overheating  occurs  will  be 
determined  and  analized. 

In  order  to  calculated  the  maximum  value  of  the  product  of  u;t,  the  electron  cyclotron 
frequency  times  the  average  collisiontime  (tc  ~  £)  in  th  constrictor  of  a  nozzle  type 
hybrid  thruster  (see  Fig. 1.1a),  one  has  to  determine  at  first  the  condition  within  the 
conctrictor  region  where  the  magnetic  field  becomes  a  maximum  (see  Fig.3.1). 

By  modelling  the  rotational  symmetric,  axial  current  density  distribution  within  the 
constrictor  through  a  paraboloid  of  grade  n  (see  Appendix  A.2)  one  obtains  for  the 
azimuthal  selfinduced  magnetic  induction  field 


B$  =  B(r)  =  /i0^  J  jrdr 

ftp  I  n +  2  r  r  2  /r\"l 
2t  rc  n  rc  L  n  +  2  \rc)  J 


(3.1) 


The  radial  dependent,  absolute  value  of  the  electron  cyclotron  frequenty  follows  therefore 
(see  also  Appendix  A.2)  by 


we 


cfl(r) 

nt. 


(3.2) 


34 


V.v 


Figure  3.1:  Illustration  of  the  Self  Magnetic  Field  within  a  Plasma  Thruster 


And  the  maximum  value  of  Be  and  uv  follows  now  by  the  requirement 

|:ft(r)*|«.(r)=  0 

r  /I  n  +  2\« 

To  V  2  Ti  +  1  / 

Hence  one  obtains  the  optimal  gyration  frequencey  by 

=  -L  as-L  ( I\"  (“HA*  =  J-fSL Lj3 

me2xrc\2/  Vn+1/  m,2rrc 


(3.3) 

(3.4) 

(3.5) 


where  J3  is  a  function  of  the  exponent  n  of  the  current  density  paraboloid  which  varies 
between  about  1.1  and  1.0  if  n  changes  from  1.5  to  100. 

Now  the  electron  collision  time  re  or  the  reciproke  value  of  the  collision  frequency  v~l 
follows  according  to  Appendix  A.l  by 
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Since  Qtt  can  be  expressed  by  the  Grosdover  crossection  [see  Appendix  A.2,  eq.  A.59 
with  2„  =|  znt  |=  1]  and  since  the  sum  with  the  terms  proportional  to  are 

negligible  with  respect  to  one,  the  collision  term,  re,  for  the  electrons  reduces  to 

1  24v/5*i 

"-^{1+144^}  ^ 

(3.7) 

By  means  of  the  definition  of  the  degree  of  ionization 

a  = - = - c - 

ft*  ft* 

(3.8) 

and  the  following  relations 

p  =  5^  n„kT  +  nekTe 

vft 

One  obtains  with  Te  =  T„  =  T  (constrictor  condition) 

(3.9) 

a  P 
'  1  +  a  kTc 

(3.10) 

and  can  replace  n,  in  the  formula  for  tc  by  p,  Tc  and  a. 

The  dimensionless  product  u,re  which  is  of  interest  here  follows  now  from  eq.  3.5,  3.7 
and  3.10  by 


12\Z2rJ3  po«o  (kTr)i  I 
W'T'  “  /„{i  +  72,r^}e3v^7  nc  rc 

__  _ dy/2vJ3 _ Pofp  l+q(U,)^/ 

/n|5.7-  10-3T,(^)i|e3\/fS'  °  p  Tc 
1.337  10 -13j3  1+0(7;)*/ 

/n|5.7  10-37;(^)^  Q  p  r' 


were  the  electron  temperature  Tt  must  be  given  in  [K],  the  pressure,  p  in  [Pa]  and  the 
quotient  arc  current,  I,  divided  by  arc  radius  rc  in  [Aj .  The  degree  of  ionisation,  a  is 
herein  a  funktion  of  Tz  and  p.  The  electron  temperature  Te  within  eq.  3.11  follows  now 
by  that  at  the  channel  radius  given  by  eq.  3.4.  The  temperature  profile  is  now  according 
to  Appendix  A.2 


Te(r)  =  T, ; 


n  +  2  fr\ 2  t-Errspfe)  ~W> 

„(„  +  4)W  l-$$<(r)  J 


(3.12) 


where  n  is  the  exponent  of  the  current  density  paraboloid  and  (( r )  ~  £(r)  =  (  is  a 
dimensionless  factor  which  accounts  for  radiation  losses  and  which  can  be  considered 
small  compared  to  one. 
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n  = 

1.5 

2 

3 

4 

10 

100 

0 

II 

w 

0.701 

0.680 

0.646 

0.619 

0.522 

0.292 

0.05 

0.690 

0.676 

0.643 

0.617 

0.521 

0.292 

0.1 

0.690 

0.671 

0.640 

0.614 

0.520 

0.291 

0.2 

0.673 

0.659 

0.632 

0.608 

0.517 

0.291 

0.4 

0.588 

0.607 

0.602 

0.587 

0.510 

0.291 

Table  3.1:  J2  according  to  eq.(3.12)  and  (3.13) 


Hence  one  obtains  for  the  electron  temperature,  Te  in  eq.  3.10 


Te  =  Mn,()Tt 


(3.13) 


where  Te  is  the  maximum  centerline  temperature  and  where  the  factor  J 2  follows  by 


Ji  = 


(n  +  nn+2\i 
n(n  +  4)  \2n  +  1/ 


(3.14) 


and  is  tabulated  in  table  3.1  for  n  between  1.5  and  100  and  (  between  0  and  0.4. 


according  to  Appendix  A.2  [eq.  A.47],  one  can  write  the  maximum  electron  temperature, 
Tt  in  the  form 

t  =  1.51  •  10Vo(;r)*  (3.15) 

where  7  •  Jo  the  product  of  two  faktors  which  varries  slowly  with  ft  or  with  (£),p,n 
and  depending  on  the  type  of  propellant  (see  Appendix  A.2  for  hydrogen  and  argon). 
The  faktor  J0  is  dimensionless  and  depends  only  on  n  and  £  [see  table  B,l]  while  the 
dimensionless  7  depends  weakly  on  p  and  Tt  and  the  type  of  propellant.  Both  faktors 
?  o  of  the  order  one.  In  eq.  (3.14)  Te  follows  in  Kelvin  if  I  and  rc  are  taken  in  Ampere 
.nd  A'eter,  respectively.  If  one  now  replaces  Tt  in  eq.  3.10  according  to  eq.  3.12  and  eq. 
3.14  one  obtains 

(>•>«) 


where 


3.75  - lQ-g(J0Ja)^  1  +  a  j  r  AH 
In  |5.7  10-3Jaft(^)^|  Q  ^  ^ 


(3.17) 


is  a  coefficientof  the  order  10-8  [$]  which  varries  slowly  with  Tt  or  £  and  with  p,n  and  ( 
depending  on  the  type  of  propellant.  The  faktor  #  is  tabulated  for  hydrogen  and  argon  in 
table  3.2  and  table  3.3,  respectively,  for  different  pressure  p  and  maximum  temperatures, 
7i,with  n  =  2,4, 10  and  (  «  0. 
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T'=\ 

2 

2.5 

3 

3.5 

4 

4.5 

5 

•lOnA'] 

P= 

n=2 

3.37 

3.46 

3.62 

3.74 

3.92 

4.06 

4.19 

[$] 

103  Pa 

4 

2.44 

2.31 

2.40 

2.51 

2.60 

2.69 

2.77 

•IO"8 

10 

2.54 

1.41 

1.40 

1.45 

1.50 

1.54 

1.59 

2 

3.59 

3.68 

3.82 

3.96 

4.09 

4.21 

10 4Pa 

4 

2.50 

2.48 

2.56 

2.56 

2.73 

2.81 

•10~8 

[$] 

10 

1.96 

1.51 

1.50 

1.54 

1.58 

1.62 

2 

4.33 

3.96 

4.01 

4.11 

4.23 

4.34 

10*Pa 

4 

3.53 

3.78 

2.72 

2.77 

2.84 

2.91 

•IO"8 

[#] 

10 

4.14 

2.11 

1.72 

1.66 

1.68 

1.71 

n=2 

5.17 

4.62 

4.55 

4.59 

4.66 

106  Pa 

4 

4.21 

3.36 

3.16 

3.14 

3.17 

■10~8 

l#l 

10 

4.51 

2.71 

2.15 

1.98 

1.94 

Table  3.2:  0  for  Hydrogen 


P= 

103  Pa 


1(P  Pa 


10s  Pa 


106  Pa 


t  = 

n=2 

4 

10 

2 

4 

10 

2 

4 

10 

2 

4 

10 


2  2.5  3  3.5  4  4.5  5 

4.45  4.37  4.51  4.57  5.12  5.82  6.58 

4.04  2.99  3.27  3.23  3.44  3.87  4.37 

7.21  2.21  1.99  2.18  2.21  2.28  2.52 

4.37  4.44  4.37  4.39  4.74  5.22 

3.37  3.07  3.19  3.11  3.22  3.49 

3.93  2.12  2.00  2.10  2.11  2.12 

5.77  4.32  4.18  4.02  3.90  4.01 

5.63  3.28  2.95  2.92  2.82  2.78 

8.92  3.31  2.11  1.89  1.88  1.87 

5.47  4.03  3.61  3.35  3.12 

5.05  3.17  2.63  2.43  2.28 

6.18  3.20  2.07  1.67  1.53 


■io,W 

•io-  ISl 


Table  3.3:  0  for  Argon 


The  degree  of  ionisation,  a  in  eq  3.17  has  been  calculated  for  one  two  and  three  times 
ionized  particles  under  the  primise  that  locally  the  Saha  equilibrium  is  valid.  At  lower 
temperatures  and  higher  pressures,  a  =  a(Tt,p)  is  fairly  amall  (<  1)  such  that  0 
becomes  quite  large  (>  8).  This  range,  however,  is  not  of  interest  here  and  would  result 
in  wrong  0—  values  since  under  those  conditions  the  second  term  in  the  denominator  of 
eq.  3.6  can  not  be  neglected  anymore.  Since  here  only  the  main  electrically  conducting 
plasma  channel  is  considered,  the  space  for  the  0-  value  in  table  3.2  and  table  3.3  at 
lower  temperatures  and  higher  pressures  are  left  blanc. 

One  should  remark  here  that  eq.  3.10  and  eq.  3.11  and  therefore  eq.3.13  and  eq.3.14 
has  been  derived  under  the  premise  (see  Appendix  A.2)  that  the  radial  heat  conduction 
of  the  electrons  are  still  not  affected  by  the  u)ere—  effects;  i.e.  one  has  assumed  within 
the  energy  equation  for  the  electrons  that  the  heat  conduction  coefficient  is  a  scalar  (see 
eq  A.39]  and  that  the  factor  with  which  the  heatconduction  coefficient,  A,  of  the 
radial  flux  should  be  multiplied  fsee  Appendix  A.l  eq.  A.31)  has  been  tacitly  set  equal 
to  one.  For  small  were-values  this  assumption  is  justified  since 


sinu>eTe 


«  1  -  g(wcT,)J  +  ... 


(3.18) 


Now,  for  <  0.25  the  error  in  the  heatconduction  coefficient,  Ae  amounts  to  <  1% 
for  u>tTc  <  0.75  the  error  amounts  to  about  10%  and  less.  As  long  as  stays  well 
below  about  0.75,  one  therefore  may  apply  eq.  3.15  to  determine  u,rt.  On  the  other 
hand  the  value  calculated  according  to  eq.  3.15  schould  never  exceed  a  critical 
value  (u VTeJcrK  <  0.75  in  order  to  avoid  ”  run  away  heating  "  within  a  discharge  channel 
radius  given  by  eq.3.4.  If  for  instance  were  would  exceed  such  a  critical  value,  (h;(rc)erit> 
the  radial  heatconductionloss  of  the  electrons  would  be  cut  off  to  a  certain  extend  and 
the  temperature  within  this  core  region  should  increase.  According  to  eq.  3.1 1  also  u>cTe 
should  increase  futher  which  again  leads  to  an  even  higher  temperature  and  so  on.  This 
run  away  heating  effect  presents  an  instability  of  the  discharge  and  eventually  causes 
the  arc  core  channel  to  assume  a  configuration  by  which  also  an  increased  heatloss  is 
possible.  Instead  of  a  straight,  rotational  symetric  channel  configuration,  the  discharge 
core  region  may  become  bent  (e.g.  helically  shaped)  thereby  experiencing  an  increased 
cooling  due  to  transverse  magnetic  pumping  effects.  The  overheating  in  the  disturbed 
core  region  is  therefore  balanced  by  an  increased  cooling.  In  reality  one  observes  an  onset 
of  voltage  fluctuations  and  an  increase  in  average  voltage  at  a  certain  critical  current 
which  could  be  explained  by  such  a  run  away  heating  effect.  The  conditions  under  which 
onset  occures  are  now  mathematically  given  by  the  inequality 


W«Te  =  0-  (— )  <  (WeT.jcr*  (3.19) 

p  Vrc/ 

where  the  critical  value  should  be  at  least  large  enough  to  have  an  effect  on 

the  heatconduction  (>  0.1)  and  smaller  than  0.75  since  only  then  the  equality  sign  in 
eq.  3.16  is  valid.  One  therefore  may  limit  (w.Te)^  by  about  the  following  range 


0,1  <  (weT,)ert(  <  0.75 


(3.20) 
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Instead  of  taking  the  pressure,  p,  as  variable  in  the  stability-  or  onset  criteria  [eq.  3.17), 
one  may  replace,  p,  by  the  mass  flow  rate  m  of  the  thruster,  provided,  that  the  discharge 
fills  the  entire  cross  section  (xr|)  of  the  constrictor.  One  now  obtains  for  the  mass  flow 
rate 

m  =  2t  /  pvrdr  (3.21) 

Jo 

where  p  is  the  local  density  and  given  by  the  relation  [see  also  eq.  3.9  with  T,  =  T\ 


p  =  p(l  +  a)£-r 

mo 


(3.22) 


with  m0  as  atomic  mass  and  where  the  axial  velocity,  v,  can  be  replaced  by  the  local 
speed  of  sound  under  isothermal  conditions  [41];  hence  it  is  according  to  eq.3.20 


v  = 


'  T=conat 

The  axial  mass  flow  rate  follows  therefore  by 

P 


(l  +  o)£-r 

m  o 


(3.23) 


m  =  2ir 


r _ 

Jo  .r±A 


-rdr 


2  !mo 
=  *rcP\lVr 


1 


kT  1  -f  ot 


(3.24) 


Herein  one  neglects  the  magnetic  pinch  effect  and  takes  T  according  to  Appendix  A.2  . 
For  n  =  4  and  £  <  0.2  one  may  very  well  approximate  eq  3.24  by 


m  =  trip, 


mo  1.15 


kT  vT+6 

where  now  d  is  an  average  degree  of  ionisation  which  may  be  taken  as 


(3.25) 


with 


a  =  a(T) 

(3.26) 

f  «  0.75 f 

(3.27) 

Within  the  constrictor  the  electron-  and  heavy  particle  temperature  can  be  considered 
equal.  Hence  the  pressure  p  in  the  stability  criteria  eq.3.19  can  be  replaced  by 


P 


m  kTe  VT  +  q 
xr;}\  m0  1.15 


25.24 


m 


t,{  i  +  a) 


Mc 


[Pa] 

[Pa] 


(3.28) 
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where  the  last  expression  follows  according  to  eq  3.14;  Mc  is  the  atomic  weight.  The 
stability  criteria  can  now  be  written  in  the  form 


-*v4< 

m 


(WeTt)„ 


(3.29) 


with 


P 


=  3.22  ■  10~3 


\JiMi  +  «) 


kg 

smFAF 


(3.30) 


a  slowly  varying  factor  in  f  or  —  and  in  p  ,n  and  £  and  which  depends  on  the  type  of 
plasma.  p  is  tabulated  for  hydrogen  and  argon  in  table  3.4  and  table  3.5,  respectively, 
for  different  temperatures  Te  and  pressures  p  with  n  =  2,4, 10  and  £  <  0.2. 


t=  1 

2 

2.5 

3 

3.5 

4 

4.5 

5 

•10*1*1 

p= 

n=2 

6.29 

6.03 

5.46 

5.27 

5.36 

5.47 

5.57 

103Pa 

4 

4.70 

4.15 

3.70 

3.64 

3.70 

3.77 

3.84 

•io-n 

I  **  u  1 

10 

5.05 

2.61 

2.22 

2.18 

2.21 

2.24 

2.28 

2 

6.66 

6.34 

5.79 

5.65 

5.71 

5.79 

10"Pa 

4 

4.82 

4.40 

4.00 

3.80 

3.97 

4.02 

TO-11 

1  k3  nl 

10 

3.90 

2.75 

2.41 

2.36 

2.38 

2.40 

2 

9.04 

7.47 

7.02 

6.48 

6.35 

6.22 

b&r] 

105  Pa 

A 

7  *Q 

*  A  A 

4.91 

4.50 

4.35 

4.34 

-io-u 

10 

9.13 

4.26 

3.18 

2.77 

2.66 

2.64 

2 

11.1 

9.03 

8.31 

7.72 

7.28 

106  Pa 

4 

9.32 

6.79 

5.96 

5.44 

5.12 

■io-u 

(  *s  J 

10 

9.41 

4.40 

4.16 

3.52 

3.22 

Table  3.4:  0‘  for  Hydrogen 
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p= 

103  Pa 


104  Pa 


P= 

KPPa 


106  Pa 


n=2 

4 

10 

2 

4 

10 

n=2 

4 

10 

2 

4 

10 


5.07  4.56  3.98  3.69  3.96  4.36  4.77 

4.76  3.23  2.97  2.71  2.71  2.77  3.30 

8.77  2.45  1.86  1.89  1.85  1.84  1.97 

4.97  4.63  3.96  3.69  3.84  4.11 

3.97  3.30  2.98  2.72  2.72  2.86 

4.79  2.33  1.92  1.90  1.84  1.80 

7.55  5.10  4.56  3.92  3.54  3.50 

7.57  4.01  3.31  2.94  2.65  2.53 

5.00  4.18  2.43  1.95  1.83  1.76 

7.64  5.17  4.35  3.73  3.23 

7.29  4.21  3.27  2.79  2.44 

10.0  4.38  2.64  1.96  1.69 


■  lO-io 


■  lO-io 


■lO-io 


•  lO-io 


Table  3.5:  /?'  for  Argon 


Constrictor  diameter 

Measurement 

m  Icrit 

u>cTe  calculated 

n=10  n=4  n=2 

0.024 

0.3  2700  2.07 

0.6  3800  1.91 

0.8  4500  1.94 

1.0  4800+  1.75 

1.5  6000+  1.74 

0.38  0.57  0.78 

0.35  0.55  0.80 

0.37  0.56  0.77 

0.31  0.50  0.70 

0.35  0.49  0.71 

0.030 

0.8  4300  1.87 

1.6  6400  1.92 

0.34  0.56  0.75 
0.34  0.54  0.78 

M 

■10~3[Jkff/s]  [A]  -109 

0.35  0.54  0.76 

Table  3.6:  were- Values  calculated  for  different  Onset  conditions  at  n=2,  4,  10  for  Argon 


Constrictor  diameter 

-i,  r  . 

m  Icr.t  —  J 

n=10  n=4  n=2 

0.024 

0.15  1400  3.07 

0.1  0.15  0.17 

M 

•10_3(Atf/,s}  [A]  -109  m}lhsA9,ilkg 

Table  3.7:  uv  revalues  calculated  for  different  Onset  conditions  at  n=2,  4,  10  for  Hydro¬ 
gen 


In  table  3.6  and  3.7  the  u> cre  -  values  for  two  MPD-thrusters  with  different  constrictor 
sizes  and  propellants  are  calculated  for  the  given  onset  conditions.  In  the  case  for  Argon 
these  values  lie  between  about  0,35  and  0,76  depending  on  the  exponent  n  for  the  current 
density  paraboloid.  Hence  it  is  shown  that  the  observed  onset  occures  at  u>cre  -  values 
which  cannot  be  neglectable  besides  one  and  that  the  critical  value,  (u>c7v)cr;<  at  which 

onset  occures  is  almost  proportional  to  ^J-  what  again,  agrees  very  well  with  the  expe¬ 
rimental  results.  The  one  hydrogen  measurements  available  with  the  24  mm  constrictor 
results  in  a  somewhat  smaller  critical  value  (weTe)CTif  than  that  calculated  for  Argon. 
This  discrepancy,  however,  could  be  explained  by  the  fact  that  the  discharge  channel 
does  not  fill  the  constrictor  cross  section.  This  is  the  more  likely  since  the  volume  gas 
flow  for  hydrogen  is  much  higher  than  that  for  argon  in  a  simmilar  current  range  and 
therefore  confirms  the  current  carrying  cross  section  for  hydrogen  within  a  smaller  area 
than  that  for  argon.  The  consequence  of  this  effect  means  that  the  taken  mass  flow  rate 
for  hydrogen  is  too  high.  Nevertheless  this  calculation  shows  that  the  u>ere-values  in  a 
discharge  column  becomes  large  enought  to  explain  the  onset  phenomenon  by  ”  run  away 
Joule  heating  "  . 
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4.  Electrode  Effects 


4.1  Cathode  Experiments 


The  scheme  and  test  assembly  of  the  cathode  erosion  experiment  are  shown  in  Fig.  4.1 
and  4.2, respectively.  The  cathode  consists  of  a  thoriated  (2%  thorium  oxide)  tungsten 
rod  of  3  mm  diameter.  The  anode,  a  metal  button  of  the  same  material,  is  mounted 
transversely  to  the  cathode  rod  in  a  distance  of  1mm  from  the  cathode  surface.  The 
electrode  system  and  the  ignitor  are  located  in  the  middle  of  a  stainless  steel  tank  with 
a  diameter  of  1  m  and  a  length  of  2  m. 

By  means  of  a  PFN  (pulse  forming  network)  -  battery  and  an  ignitor  the  cathode  is 
charged  by  a  fairly  rectangular  current  pulse  of  about  1400  A  for  «  2  ms.  A  typical 
current  time  pulse  is  plotted  in  Fig.  4.3.  One  pulse  transfers  an  electric  charge  of  about 
2,6  As  (1  As  =  1  Ampere  second  =  1  Coul).  Small  variations  in  the  pulse  profile  were 
averaged  over  50  shots  per  measuring  point  and  the  small  variation  in  the  electrical 
charge  per  pulse  have  been  accounted  for. 

The  repetition  rate  of  the  fully  automatically  working  test  facility  can  be  varied  between 
several  seconds  and  several  minutes.  In  the  cold  cathode  experiment  the  repetition  rate 
is  taken  between  1  and  2  minutes  in  order  to  allow  the  cathode  to  cool  down  between 
the  shots. 


4.2  Measurements 


The  weight  loss  of  the  cathode  sample  has  been  measured  after  50  shots  by  means  of  a 
very  accurate  scale  and  the  erosion  rate  given  by  the  wheight  loss  per  transferred  electric 
charge  is  plotted  in  Fig.  4.4  a,  b,  c,  d,  e  and  f  for  ambient  pressures  of  10,  100,  1000, 
10000  and  100  000  Pa,  respectively.  These  measurements  indicate  a  strong  increase  of 
the  erosion  rate  with  an  increasing  number  of  shots  or  with  an  increasing  accumulated 
electric  charge  in  the  beginning  and  after  about  500  shots  (each  measuring  point  repres¬ 
ents  =»  50  shots  )  which  corresponds  to  about  1300  As  the  erosion  rate  levels  off  and 
reaches  a  constant  value.  This  fact  may  be  caused  by  roughening  of  the  original,  polished 
surface  due  to  spot  and  crater  formation  and/or  due  to  depletion  of  the  thorium  oxide  at 
the  main  impact  area  which  directly  faces  the  anode.  The  measurements  at  1  bar  (Fig. 
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Figure  4.1:  Scheme  of  Cathode  Testrig 


for  the  heating  current ) 


‘■iir.vs.-jVgeV&t 


t 

i 

(. 


Figure  4.3:  Typical  Current  vers.  Time  Curve 

4.4f)  scatter  somewhat  since  the  weight  loss  in  that  case  is  fairly  small  («  0.05  mg)  and 
therefore  difficult  to  measure  correctly.  In  Fig.  4.5  these  results  are  again  summerized 
in  one  plott  in  order  to  show  the  influence  of  the  pressure  on  the  erosion  rate. 

If  one  plots  the  asymptotic  value  of  erosion  rate  as  funktion  of  pressure  as  shown  in  Fig. 
4.6,  there  is  a  distinct  threshold  between  an  ambient  nitrogen-  pressure  of  100  and  1000 
Pa.  Below  about  100  Pa  the  erosion  rate  seems  to  stay  fairly  constant  at  27  / ig/Coul , 
above  100  Pa  the  erosion  rate  drops  threshold  like  with  increasing  pressure  to.  about  9 
gg/Coul  at  1000  Pa  and  about  4  gg/Coul  at  10  000  Pa  (=0.1  bar). 

Erosion  experiments  with  heated  cathode  samples  showed  only  limited  results.  The 
mass  loss  (outgasing  and  evaporation)  due  to  auxiliar  heating  alone  are  already  fairly 
high  and  the  specific  loss  amounts  to  about  0.1  gg/Ws  as  shown  in  Fig.  4.7.  Interesting 
to  note  here  is  the  fact  that  within  a  nitrogenatmosphere  the  mass  loss  increases  rapidly 
when  the  heating  power  exceeds  about  450  W  what  corresponds  to  a  sample  surface  tem¬ 
perature  of  about  2000  K.  The  calibration  curve  [see  Fig.  4.8]  between  heating  power  and 
surface  temperature  are  based  on  pyrometer  measurements  taking  an  emission  coefficient 
for  thoriated  tungsten  of  c  =  0.45.The  strong  increase  of  massloss  due  to  heating  alone 
in  a  nitrogen  atmosphere  above  sample  temperature  of  2000  K  indicates  that  chemical 
reaction  of  nitrogen  and  the  thoriated  tungsten  surface  will  enhance  the  mass  loss  of 
the  sample.  These  mass  loss  due  to  pure  heating  could  not  be  avoided  even  after  longer 
heating  periods  and  outgasing.  They  hampered  the  proper  erosion  measurements  caused 
by  the  arc  pulses.  Nevertheless,  the  erosion  rates  on  the  cathode  samples  (  thoriated 
tungsten)  due  to  arcing  alone  has  been  determined  by  measuring  the  mass  losses  with 
and  without  arc  pulses.  These  cathode  losses  show  that  the  erosion  rates  are  about  one 
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Figure  4.4:  Measured  Erosionrates  in  [^]  as  Funktion  of  the  Accumulated  Total  Electric 
Charge  in  [As]  for  Different  Ambient  Pressures.  Each  Measuring  Point  is  the  Result  of 
50  Current  Pulses.(Cathode:  Thoriated  Tungsten  at  Room  Temperature) 
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Figure  4.5:  Erosion  Rate  as  Funktion  of  the  Accumulated  Charge  (  Number  of  Current 
Pulses) 


PH£S5Uk£  iPa] 

Figure  4.6:  Asymptotic  Erosion  Rate  as  Funktion  of  the  Ambient  Pressure  of  a  Tungsten 
Cathode  in  a  Nitrogen  Atmosphere  Charged  by  a  Rectangular  Pulses  of  1400  A  for  2  ms 
(Pulsfrequency  =1-2  min~l) 
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Figure  4.7:  Massloss  Due  to  Pure  Heating  in  a  Nitrogen  and  Argon  Atmosphere  (»  < 
lOOPa) 


order  higher  at  heated  cathodes  than  at  cold  ones.  In  table  4.2  the  erosion  rates  with 
increasing  number  of  discharges  are  tabulated. 


Reliable  measurements  beyond  about  2000  K  at  the  cathode  sample  could  not  be  ob¬ 
tained  with  the  present  testrig  since  fusing  between  sample  and  fasteners  occured.  At 
these  higher  temperatures  the  pressures  were  also  limited  to  about  100  Pa  and  less  since 
convection  cooling  became  too  strong  and  also  did  not  allow  to  achieve  higher  tempera¬ 
tures  than  about  2000  K. 


50 


Number  of  Discharges 

Accumulated  Electric  Charge  in  As 

50 

135.35 

31.4 

100 

270.06 

186.0 

150 

405.89 

292.0 

200 

541.03 

226.0 

250 

678.83 

223.5 

300 

815.12 

231.3 

350 

949.52 

220.9 

400 

1082.92 

243.7 

450 

1219.24 

250.5 

500 

1354.68 

229.8 

average 

2.71 

scattering 

2.14  % 

Average  Asymptotic 
Erosionrate 

245.2 

Table  4.1:  Erosion  rate  caused  by  arcing  only  at  the  surface  of  a  heated  tungsten  (2% 
thoriated)  cathode  in  a  Nt-  atmosphere  (<  lOOPo)  surface  temperature:  «  1500 A' 


4.3  Conclusion  and  Explanation  (qualitative) 


The  falling  characteristic  of  the  erosion  rate  with  increasing  pressure  (  see  Fig.  4.6  )  has 
been  observed  before  [40,  41,  42,  43]  and  reported  on  .  Here  a  qualitative  explanation 
based  on  the  previously  developed  spot  or  micro  jet  theory  [18,  38,  44,  46,  47,  48,  49]  is 
presented.  Based  on  that  theory  the  main  electron  emitting  cathode  area  consists  of  one 
or  more  hot  spots  of  high  current-density  (j  as  lO'M/m2  and  higher)  which  cause  locally 
an  extremely  high  heat  load  on  the  cathode  surface  leading  to  strong  local  evaporation 
and  micro  craters  from  which  electrically  conducting  plasmajets  propagate.  Within  the 
spotcraters  the  vapor  plasma  pressure,  p,,  is  quite  high  and  can  reach  several  tens  to 
even  100  bar  which  depends  on  the  electrode  material  while  the  ambient  pressure, p<», 
is  usually  much  smaller.  Now,  the  spot  theory  tells  us  that  a  spot-  discharge  remains 
stable  i.e.  sticks  to  the  same  crater  like  emission  site  as  long  as  the  requirement 

P.  -  Po  >  3.3  g/.j,  (4.1) 

is  fulfilled  i.e.  the  pressure  at  the  inner  crater  surface  p.  minus  the  pressure  at  the  outer 
edge  of  the  current  carrying  plasma  jet,  po,  (  see  Fig.  4.9  )  is  larger  than  about  3.3-times 
the  average  pinchpressure  of  the  electrically  conducting  channel  of  the  plasmajet.  If 
the  requirement  (4.1)  does  not  hold  anymore  the  main  current  carrying  channel  of  the 
plasmajet  will  bend  and  attaches  the  cathode  surface  somewhere  in  the  neighborhood  of 
the  original  spot  thereby  creating  a  new  spot  site.  Therefore  the  spot  current,  /,  cannot 

p,"2sr"o  Jet 


Figure  4.9:  Scheme  of  a  Spot  Discharge 
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Figure  4.10:  Illustration  of  the  Plasmaflow  around  the  Crater  Rim  with  Prandtl- Meyer 
Deflection  Angle, 6 

exceed  a  value  as  given  by  eq.  (4.1). 

In  addition  the  spot  theory  shows  that  the  minimum  erosion  rate,  eu,„  in  [£j],  due  to 
evaporation  is  independent  of  the  total  current,  /(=  /,),  or  the  number  of  spots, 

but  depends  on  the  (average)  spot  current  /„  the  molecular  or  atomic  weight,  M  of  the 
evaporated  cathode  material  and  the  surface  temperature  within  the  active  crater,  T,. 
It  is 

CMm  =  c^r/.  <e  (4.2) 

where  e  is  the  real  erosion  rate  and  |c  ss  10_9[^rr]j  ’s  a  ^actor  which  depends  in  a 
minor  way  on  the  geometry  of  the  crater  and  to  some  extent  on  an  applied  magnetic 
field.  Here  in  our  consideration  in  which  the  directional  motion  is  omitted,  one  may 
take  c  as  a  constant.  Based  on  the  relations  (4.1)  and  (4.2)  one  can  now  explain  the 
threshold  like  falling  characteristic  of  the  erosion  versus  pressure  curve  by  means  of  the 
Prandtl-Meyer  deflection  [50]  as  shown  in  Fig.4.9  and  4.10. 

According  to  eq.  (4.2)  the  minimum  erosionrate  and  consequently  the  real  erosionrate 
decreases  the  smaller  the  average  spot  current,  /„  at  which  the  spot  discharge  becomes 
unstable. 

Now,  in  order  to  achieve  such  an  unstable  condition  (see  eq.(4.1)  one  has  to  make  at 
a  given  spot  current  density,  j,  the  pressure  difference  (p,  -  po)  smaller,  or  at  a  given 
pressure  difference  the  current  density  larger.  The  average  current  density,  j,  however,  is 
strongly  coupled  to  the  spot  pressure,  p,,  and  both  are  constant  or  simultaneously  either 
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increase  or  decrease.  Both  quantities  j  and  p,  cannot  be  changed  directly.  Therefore, 
in  order  to  lower  the  spot  current  and  the  erosionrate  one  has  to  increase  po  or  better 
increase  the  ambient  pressure,  p« .  This  fact  agrees  qualitatively  with  the  experimental 
finding. 

The  threshold  like  increase  of  the  erosionrate  by  lowering  the  pressure  from  about  100 
Pa  to  1000  Pa  (see  Fig.  4.6)  needs,  however,  a  more  detailed  discussion  of  the  shape  of 
the  plasmajet. 

As  long  as  the  Prandtl- Meyer  deflection  6  is  smaller  than  a  critical  angle,  8^,,  (see  Fig. 
4.10,  where  8crit  spans  between  the  slope  of  the  inner  crater  surface  and  the  proper 
cathode  surface)  the  edge  of  the  plasmajet  does  not  touch  the  cathode  surface  outside 
the  crater  and  the  pressure,  po  can  be  considered  equal  to  the  ambient  pressure,  Poo-  If, 
however,  the  Prandtl- Meyer  deflection  exceeds  6„u,  the  plasmajet  will  hit  the  cathode 
surface  and  thereby  evaporate  cathode  material,  (po  >  poo)  and  eventually  increase  the 
crater  size. 

The  Prandtl- Meyer  deflection  angle,  8  increases  now  with  decreasing  pressure  ratio,  ®-, 
according  to  the  following  relations  [50] 

6  =  tfi  -ip  (4.3) 


where  ip  and  ip  follows  from 


a  + 1  po 


(^)^-l 


Un(vyO),^ 


(4.4) 


(4.5) 


k  =  is  the  ratio  of  the  specific  heats,  with  the  range  of  the  pressure  ratio:  [o  <  ^  <  {^ff}'r^T 

For  (=0.56  for  k=1.2  and  0.49  for  k  =  1.66)  the  expansionflow  reaches 

Ma  =  1;  for  &  <  (^~y)“^T  the  expansionflow  is  supersonic.  The  deflectionangle  8  is 
tabulated  for  various  values  m  and  k  in  Table  4.2. 

P» 

For  an  ambient  pressure  p«,  of  about  I6ar(=  ltPjPa)  one  can  expect  the  pressure  ratio 
“■  in  the  range:  =  ?  =  10_1  -  10-2;  i.e.  the  Prandtl-Meyer  deflection  can  still  be 

considered  small  such  that  the  plasmajet  does  not  hit  the  cathode  surface  outside  the 
crater  rim.  If,  however,  reaches  values  of  about  10-3  —  10~4,  there  is  a  great  possibi¬ 
lity  (even  for  k  =  1.66)  that  the  plasma  comes  close  to  or  even  hits  the  cathode  surface 
and  therefore  evaporates  additional  cathode  material  (po  >  pm).  By  further  lowering  of 
the  ambient  pressure,  p this  Unfitted  area  will  also  emit  electrons  and  must  be  counted 
to  the  spot  crater.  Hence  the  spot  size  of  the  electrically  active  crater  increases  if  one 
lowers  the  ambient  pressure.  This  fact  has  been  indeed  observed  before  [51].  By  means 
of  the  well  known  phenomenon  of  the  Prandtl-Meyer  expansionflow  applied  to  the  flow 
of  cathode  spot,  one  therefore  can  explain  the  almost  thresholdlike  increase  of  the  ero¬ 
sionrate  by  lowering  the  ambient  pressure  below  ~1000  Pa.  The  observed  erosionrate 
versus  pressure  curve  as  shown  in  Fig.4.6  ,  however,  is  specific  for  the  given  experiment 
and  cannot  be  generalized.  At  more  powerful!  current  pulses  and/or  longer  lasting  pul¬ 
ses  than  in  the  experiment  presented  here  and  especially  on  ductile  cathode  materials 
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Table  4.2:  Prandel-  Mayer  Deflection  6  for  various  pressure  ratios  21  and  k-  values 


like  copper  and  aluminium  etc.,  the  transition  from  higher  erosionrates  at  low  ambient 
pressures  (<  lOPa)  to  lower  erosion  rates  at  higher  ambient  pressures  (>  103  Pa)  is  less 
steep  and  moreover  the  erosionrate  can  increase  dramatically  at  higherpressures  too. 
This  latter  behavior  stems  from  the  fact  that  at  higher  pressures  (>  104  Pa)  the  spot 
activity  is  located  in  a  smaller  area  of  the  cathode  surface  and  that  not  only  the  ac¬ 
tive  spot  but  the  entire  area  becomes  molten  and  erodes  very  heavily  by  droplet  jection 
and/or  splashing  [45]. 
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A.  Appendices 
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A.l  Influence  of  the  Magnetic  Field  on  the  Heat- 
flux  within  an  Electrically  Charged  Plasma- 
component 


It  has  been  shown  previously  1  that  in  a  relaxation  free  plasma  (distribution  functions 
of  different  gascomponents  adjust  themselves  immediately  to  a  given  special  condition) 
the  distribution  function  fj®  due  to  gradients  becomes  pertubed  and  can  be  written  in 
first  approximation  by 

fM  ~  fl°)  _  .  v/«°>  (A.l) 

where  //°*  is  the  Maxwell  Boltzmann  distribution 


and  A£,  represent  a  vector  elemement  with  the  dimension  of  a  length  the  absolute  value 
of  which  corresponds  to  a  mean  free  path  i.  e.  the  probalitity  length  along  which  a  test 
particle  sustaints  and  transports  its  velocity  V<  or  its  distribution  properties.  In  that 
derivation  one  has  assumed  so  far  that  no  external  field  effects  are  present.  If  one  now 
allows  that  a  magnetic  induction  field  B  is  present,  one  has  to  account  for  the  fact  that 
the  ”  statistical  motion”  of  aZ,e  -  charged  particle  proceed  not  anymore  along  a  straight 
line  but  rather  along  a  circular  or  helical  line  caused  by  the  magnetic  force  Z,e(Vi  x  B) 
2.  In  the  following  one  now  consideres  a  magnetic  field  B  perpendicular  to  the  gradient 
of  the  distribution  function  /,■ 0) . 

In  the  presents  of  a  magnetic  field  it  is  therefore  convenient  to  express  the  length  vecto- 
relement  by  the  sum  of  two  vectors  :  A^u  which  is  paralell  and  A£, j.  which  is  perpendi¬ 
cular  with  respect  to  B.  Hence 

A&  =  A|*j|  +  a£x  (A.3) 


Within  a  magnetic  induction  field  B  a  test  particle  with  the  electric  charge  Z,e  and  the 
peculiar  velocity  Vi  —  V^|  +  Vfj.  (  ||  parallel,  ±  perpenticular  with  respect  to  B  )  moves 
along  a  circular  or  helical  line  with  a  Lamorradius  given  by 


,  mj  [%  x  B\ 

rL= - — 

Zie  B 2 


I  n  1=  rL  = 


m,  1  | 

Zit  B 


(A.4) 


the  statistical  pathlength  is  now  given  by 


I  Af  |=  As 


^(Asx)2  +  (As,,)2 


(A.5) 


'Interim  Scientific  Report  "Basic  Processes  of  Plosmapropulsion”  Grant  APSOR  86-0337,  covering 
period  l.Aug  1988  -  31  July  1989 

Jone  may  here  approximate  the  test  particle  velocity  vj  by  its  peculiar  velocity  \7,  as  long  as  the 
average  thermal  speed  is  much  larger  than  the  relative  drift  velocity  between  the  colliding  particles 
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where 


and 


Asx  =  J-IKxl 

A*!!  =  h  I  Kill 


\ZmvkTu 
m,  +  m„ 


(A-6) 

(A.7) 


is  the  average  collision  frequency  of  the  test  particle »,  and  Qtu  the  effectiv  collision  cross 
section  between  particle  i  and  v.  According  to  eq.  (A.4),  it  follows  from  Fig.(A.l) 


Asx  =  2  rup 


(A.8) 


and  for  the  secant,  A  a  =|  A^j.  |  one  obtains 


Aa  =  2rL  j  sin  v?  | 


(A-9) 


Based  on  eqs.  (A.8),  (A.6)  and  (A.4)  it  follows 

2  =  ^ax  _  Z'eB 
rL  rmvi 

where  cj,  is  the  cyclotron  frequency  of  the  particle 


=  U),T, 

i 


(A. 10) 


and  r,  is  the  average  collision  time 


oi,-  = 


ZitB 

m, 


(A- 11) 

(A.12) 


The  length  vector  element,  A&,  in  the  perturbation  term  —  A £•  V/,-0)  can  now  be  written 
by 


1 


Vi 


Kh  +  ACx 


(A.13) 


where  the  length  vector  element  normal  to  B  (see  Fig.  (A.l)  follows  by 

A&l  =  X-  2  fL  (A- 14) 

Now  it  is 

‘A15> 

where  the  sign  of  sin  =  cos  <p  alternates  between  plus  and  minus  depending  on 

lies  within  a  range  n*  <  ip  <  (n  +  l)x  with  even  (n  =  0,2,4...)  or  uneven  (n  =  1,3,5...) 
numbers  for  n  ,  respectively.  Together  with  eq.(A.4)  and  aq.(A.9),  one  obtains 

(AJ,> 
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In  order  to  eliminate  A£tl  on  the  right  side  of  this  equation,  one  takes  the  vector  product 
of  the  equation  and  obtains 


(A&x  xB}  = 


|  siny?  I  Z,e 


(A.17) 


Indrocuced  into  the  equation  of  A£x>  one  obtains 

-  nSi  {-^?  •  BAfa  -2,^1  *x  X  B) 

|sin<£|fl[  |  sin  y?  |  Z,e  J  ZitB 2 

and  solved  for  A&±,  one  finally  ends  up  with 

A&.  =  ±2  |  sin  <fi  [cos  ^Vix- 2 sin2 ifi  ^j[V,±  x  B] 

=  sin  2v>^Kx  -  2(1  -  cos  2<p)^fej[Vix  x  B) 

or  together  with  eq.(A.lO) 

1  _  1  rsinwnp  .  i^S2«ii[Kx  *  B]} 

V  \  LJiTi  LJiTi  B  J 


(A.18) 


(A.19) 


(A.20) 
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Figure  A.l:  Illustration  of  the  Particle  Motion  within  a  Magnetic  Induction  Field,  B, 
normal  to  the  Gradient,  V/,  :  (a)  towards  V/,;  (b)  Vix  opposite  to  V/;  (B  normal 
to  plane). 
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The  entire  vector  element  is  now  given  according  to  eq.(A.20)  by 

1  .7  ,  sinu.r,  1  1-cosw.n  1  ,  - 

Mi  =  -K||  + - vi± - —  Vu.  X  B\ 

Vi  "  u\r<  Vi  UiTi  B  Vi 


(A.21) 


The  perturbation  term  in  eq.(A.l)  follows  therfore  for  the  restricted  case  in  which  the 
magnetic  induction  field  B  is  perpendicular  to  the  gradient,  V fj°\  by 


=  - 


tin  Wi  r,  1  1/  _  1—  ccww, 

U»,T, 

an  u>.  r.  1  »/  1 — coa  u/, t, 

t.l-T.  IX-  LIT. 


X  B]}  •  V/?°> 
^-[K  x  B]}  •  V/f> 


(A.22) 


Both  expressions  are  identical  since  the  directional  derivations  of  and  K  V...  and 

\V,±.  x  B]  ■  V...  and  [V*  x  B]  •  V...  give  the  same  results  respectively  (  V/  is  perpendicular 
to  B  ). 

Since  for  /,  the  following  normalization  requirement  must  be  fullfilled 


A/o,  =  /  fidVi  =  rn  (A.23) 

one  has  still  to  normalize  the  perturbed  distribution  function  /*p).  For  //p)  this  zero 
moment  becomes  not  exactly  n<  (see  the  following  side  calculation) 

MLp)  =  U\T)dVi 

—  -  f  ■  «■"(“. T.)  1  T7.r  i  l-to«:i(u,r,l  1 

-  n‘V  w.r,  w,r,  -J 

* 

with  ^ 2*77 
Now  it  is 

=  A  /  J /-+~  / G  [({.ft  -  +  (^Bv  -  evflx)£*x]  + 

+  4*  (syA*  _  -  ZxBx)^V~iy  +  (£rBy  -  + 

+  6  +  (£xfix  -  +  (exfly  -  eyflx)£»J  e-X-^dt^dZ, 

=  ^xv).«i 
=  ^[Vxuil-fl 


and  theu  fore 


coa\u>t  r, 
u»,r. 


JgtVx^B) 


Now,  the  proper  normalized  distributionfunction  /,  under  perturbed  conditions  with 
temperature  and  density  gradients  follows  according  to  this  result  from  eq.  (A.l)  by 

yf  >  -  gtot-LK-  ■  y/f >  +  x  fl]  •  V/<0) 


It  may  be  of  interest  to  note  here  that  for 

IustsI  =  I— B  <  1 

I  miVi 

the  distribution  function  /,  becomes  effectively  independend  of  B  and  assumes  the  known 
form  as  calculated  before  under  the  assumption  of  no  external  field  3 
The  heatfluxvector  due  to  conduction  follows  now  by 


£  =  /  ynWVM 


(A.25) 


By  replacing  /;  according  to  eq.(A.24)  and  summation  over  the  velocity-space,  one  even¬ 
tually  obtains  the  heat, flux  vector  due  to  conduction  within  a  magnetic  field  normal  to 
V/i  (i.e.  normal  to  V7),  Vn,)  by 

3  _  5n,(kT.)*  1 _ 

•  [^Vn,  +  £VT,|  -  i=am±  +  2^]  j 


«*  -  “ 
T‘  "  *  " 


(A-27) 


Noteworthy  is  the  fact  that  for  |  «,r,  |>  0  the  heatflux  vector  has  a  component  not  only 
parallel  to  the  temperature  and  density  gradients,  but  also  in  the  direction  normal  to 
B  and  Vn,,  VI’,.  For  increasing  |  w,r,  |  values  the  heat  flux  in  the  dirction  of  Vn„ 
VT,  dicreases  and  the  flux  normal  to  B  and  the  gradients  increases  at  first  reaches  a 
maximum  and  then  decreases.  For  |  u),t,  |=  t,2jt,3t...,  q,  diminishes.  If  one  takes 
Vn,  =  | ^VT,,  and  since  for  most  cases 


— Vn,  <  1  and  —  fv  x  e,l  <  1 
Vi  v,  1  > 


(A. 28) 


one  can  write  the  heat  flux  vector  due  to  conduction  in  good  approximation  by 

=  -KVT, 


(A. 29) 


where  now  the  heatconductioncoefficient  becomes  a  tensor  and  follows  by 


1-COIMf.  Bn 

WiT,  7 

1  —  coiu/iT,  Bm 

utiTi  B 


1  —  CO>U»,Tj  Bm 
WTi 

l-CO«U HTj  Bx 

UtiTi  B 


1— co«  U>.r.  Ba 

UfiTi  B 
l-co>u>.r.  B* 
B 


(A.30) 


3B«ic  Processes  of  Plasma  propulsion,  Interim  Scient.  Report,  Grant  AFOSR  86-0337,  covering 
period  Aug  88  -  July  89 
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with 

For 


the  heatconduction  coefficinet  becomes  a  scalar  and  independent  of  the  magnetic  field, 
B;  it  then  assumes  the  form  of  the  classical  conduction  coefficient  as  derived  elsewhere  [1] 
what  has  to  be  expected.  For  noticable  values  |  u>,r,  |>  0.1,  however,  the  heatconduction 
coefficient  is  a  tensor  the  elements  of  which  become  periodically  zero  for  lo,t,  =  ni r  with 
n  =  1,2,3... ;  and  their  relative  maximas  inbetween  become  smaller  the  larger  the  |  w.r,  |- 
values. 
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A. 2  Determination  of  the  Temperatureprofile  wi¬ 
thin  the  Constrictor 


In  order  to  calculate  the  temperature  profile  within  the  current  carrying  plasmachannel 
of  the  constrictor  (see  Fig.  A.l)  one  solves  within  the  range  0  <  r  <  rc  the  energy 
equation  of  the  electrons  under  the  assumption  that  (a)  the  axial  temperature  gradient 
can  be  neglect  besides  the  radial  one,  (b)  the  temperature  difference  between  electrons 
and  heavy  particles  like  atoms  and  ions  is  negligible  (T,  =  T)  and  (c)  the  current  density 
can  be  modelled  by  a  paraboloid  of  grade  n  according  to 


-Ha'i 

(A.32) 

where 

,  n  +  2  / 

]  = - ~ 

n  irr* 

(A -33) 

results  from  the  requirement 

.  .  .  ,  on- 

/  =  2tt  /  jrdr  =  irrt  j 

Jo  n  -f  l 

(A. 34) 

Hence  the  energy  equation  for  the  electrons  reduces  to 

J-  ==  (rA*^r)  +  RL 

a  r  dr  \  dr  J 

(A.35) 

Where  RL  stands  for  reaction  -  and  radiation  losses  of  the  electrons.  After  integration 
and  rearanging  one  obtains 

-  2 xr\'~  =  2*  ( 2{ 3-  -  RL}rdr  (A.36) 

dr  J o  <7 

Now,  within  the  arc  column  the  axial  voltage  drop  jf  or  the  electric  field  E,  can  be 
considered  independent  of  r  and  since  within  the  throat  j  =  j2  =  oEz,  it  follows 


j(r) 

<r(r) 


J 

—  =  const 

(7 


According  to  eq.(A.36)  together  with  (A. 32)  and  (A. 37)  one  obtains 


(A. 37) 


(A.38) 


The  heatconduction  coefficient,  Ae  of  the  electrons  follows  now  under  the  premise  that 
the  radial  heatflux  is  not  strongly  affected  by  u>ere  -  effects  (a vr,  <  0.1)  as 


Ae  =  Ca  •  T) 


(A.39) 
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and  the  electric  conductivity  as 


a  =  ca  ■  7V 


i 


1M 


(A. 40) 


where  c\  and  c„  are  coefficients  which  are  only  weakely  temperature  -  and  pressure 
dependent,  at  least  as  long  the  degree  of  ionisation  is  larger  than  a  few  percent. 


cA  cs 


In 


2.88  10-10  f  Te  dnA  r  W  1 
^5.7  •  10-s  J  j  l  \-mKl J 


3.88-10 


r3 


In 


(57,o'ij) 


A 

Vm/d 


d. 


(A.41) 


(A.42) 


If  one  now  replaces  in  eq.(A.38)  Ae  and  <7  according  to  eq.(A.33)  and  eq.(A.40),  resepec- 
tibely,  and  j  according  to  eq.(A.33),  one  obtains  after  integration 

f)  -  T$  =  — 


8?r2 


(A«> 

1  e  C 

«»-*};{}■ 
o  Vo  • 


where 

'  '  RL  .  ,  , 

— —rdr  dr 

.  r  o 

o  \o  J  ' 

depends  on  the  radiationloss  RL  and  can  be  considered  small  compared  with  one.  One 
can  show  that  for  a  hydrogen  plasma  at  a  pressure  below  about  5  bar  this  assumption 
is  justified,  ca  is  an  average  value  of  c\  which  is  defined  by 


-  _  T&cxridT' 
2 (f$  -  T?) 


(A.45) 


with  Te  ~  6000 A'  One  may  now  assume  that  the  temperature  Te(rc)  ~  6000A'  at  the 
edge  of  the  arc  core  is  quite  smaller  than  the  maximum  temperature  Te  in  the  center 

of  the  core,  i.e.  Te  <.  fj.  Hence  from  the  above  equation  the  maximum  temperature 
follows  with  r  =  rc  by 


8tt2  ca  Ce 


(n  +  2? 
n(n  +  4) 


(A.46) 


In  order  to  get  the  temperature  profile,  one  now  replaces  in  eq.(A.43)  the  following 
expression  according  to  eq.(A.46) 


7  /A!  _  T* _ 

Sx’cac  \rj  ~  (l  +  $  {I  - 
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One  obtains 


fi  _  T?  =  fj 


1  - 


{n±7y 


1  - 


and  solved  for  Te  yields 


r  / 1  (»  +  2)2  /r\2  f  -  pr^p-(^)  -£(r)~l 

*1  «(«  +  4)^e/  1  -  §±$«r„)  J 


(A. 47) 


The  dimensionless  temperature  profile,  £■  (—)■,  is  now  plotted  in  Fig.(A.2)  for  different 
current  density  exponents,  n.  As  can  be  shown,  the  f  -  dependency  can  be  neglected 
within  this  model  approach  as  long  as  £(r„)  ~  £  <  0.2. 

In  order  to  determine  the  maximum  electrontemperature,  Te,  one  starts  from  eq.(A.46) 
replaces  c\  according  to  eq.(A.43)  and  c„  according  to  (A.40),  respectively,  and  obtains 


t=  1.513  •  102  ■  JqI 
where  the  factor  J0  is  now  given  by 


(B  w 


Jo 


1 


(«  +  2)a 

n(n  +  4) 


«re) 


(A.48) 


(A. 49) 


and  tabulated  for  different  n  and  £  (see  Table  A.l),  and  7  is  a  dimensionless  factor  of 
the  order  one  which  follows  for  singly  ionized  gases  by 


n= 

1,5 

2 

3 

4 

10 

100 

0 

II 

1,297 

1,246 

1,185 

1,149 

1,070 

1,008 

0,1 

1,256 

1,211 

1,155 

1,122 

1,047 

0,987 

0,15 

1,233 

1,191 

1,139 

1,107 

1,034 

0,976 

0,2 

1,209 

1,171 

1,122 

1,092 

1,021 

0,964 

Table  A.l:  J0(n,  £) 


with 


(A. 50) 


(A.51) 
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and  Te  ~  6000A'.  The  numerical  factor  in  eq.  (A.48)  has  the  dimension  [A'ms/Ai], 

If  one  now  assumes  local  reaction  equilibrium  according  to  ”  Saha-Eggert”, 4  one  obtains 
the  degree  of  ionization  for  a  singly  ionized  plasma  by 

1  + - — E - \  (A. 52) 

c\T%exp{-alT)) 

where 


a  = 


_ Zvnu  f 

■it  Kv  1 


Cl 


C2 


Hi  ■  2  (27rme)*  | 
u0  h 3 
Xoi 
k 


(A.53) 

(A.54) 


The  expression  js  the  quotient  of  the  partition-  functions  and  of  the  order  one,  mc 
is  the  electron  mass,  k  is  the  Boltzmannconstant  and  h  is  Planck  constant;  xoi  is  the 
ionization  energy.  With  the  knowledge  of  a  and  the  premise  that  within  the  constrictor 
Te  =  T,  one  can  replaced  pe  and  in  eq.(A.51)  by 


pe  =  nckTc  =  -4-p  (A. 55) 

1  +  a 

and  for  a  singly  ionized  plasma  the  second  expression  by 


Tcdnc 

f. 

1., 

,  r5 

.  ^2 1 1 

nt  dTe 

—l1 

~  2  1 

-  °  [2 

+  7;!} 

n= 

1,5 

2 

3 

4 

10 

100 

4=0 

0,736 

0,743 

0.751 

0.757 

0.770 

0.778 

0,05 

0,733 

0,740 

0,750 

0,756 

0,770 

0,778 

0,1 

0,728 

0,737 

0,748 

0,754 

0,69 

0,778 

0,2 

0,716 

0,728 

0,742 

0,751 

0,768 

0,778 

0,4 

- 

0,693 

0,724 

0,738 

0,765 

0,777 

0,6 

- 

- 

0,655 

0,704 

0,757 

0,777 

Table  A.2:  ./,(«,  0 


The  integral  7*  can  now  be  calculated  for  a  hydrogen  plasma  as  function  of  Tt  and 
p(=  const)  and  finally  7,  the  expression  of  eq.(A.50).  The  somewhat  lengthy  but 
straight  foreward  numerical  procedure  leads  to  a  function  7  =  7 (p,Te)  which  is  plot¬ 
ted  in  Fig.(A.3)  for  different  pressures  in  a  temperatu rerange  between  about  10000A' 
and  50000A-.  Within  this  range  the  dimensionsless  factor  7  can  be  well  approximated  ( 
with  an  error  of  less  than  3%  )  by  the  relation 

7  =  9.35  ■  1  0~3[fe  -  2230]*  -  4.34  •  10_2/n  (A.57) 

where  fe  is  taken  in  Kelvin  and  p  in  Pascal.  Knowing  7  as  function  of  Te  and  p,  one  now 
can  determine  from  cq.(A.48)  the  maximum  electrontemperature,  fe,  as  function  of  ~, 
p,  n  and  ?(rt).  For  a  current  density  profile  given  by  a  paraboloid  with  an  exponent  of 

4Un«old,  ”  Physik  der  Stern  atmosphere  ”  ,  Springer  Book  [25] 
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n  —  4  and  a  loss  coefficient  of  £  ~  0,  the  maximum  electron  temperature,  te  has  been 
calculated  according  to  eqs.(A.48),  (A.51),  (A. 52)  and  (A.57)  as  function  of  T  and  p. 
The  results  are  plotted  in  Fig.(A.4). 


For  a  plasma  with  one  and  more  times  ionized  particles  on  can  now  determine  7  in  an 
analog  way  as  it  is  done  for  hydrogen  in  eq.(A.50).  One  obtains  now 

7  =  7*(/n 


5.7  •  10-3^f 

Pil 


i  +  E- 


q7i~z“ 


' 


(A. 58) 


where  7’  is  again  given  by  eq.(A.51).  The  ratios  of  the  collision  cross  sections  (with 
v  >  1) 


(  ^ 

'"1 

rifl44r2  £3*3f-  l 

L1  +  n«e«Zu(Z„+l)  J 

VokfJ 

>  647T 

(A. 59) 


are  mainly  functions  of  Zl  and  in  addition  to  a  minor  extend  a  function  of  the  quantities 
under  the  In  -  termes.  For  an  Argon  plasma  allowing  singly  -  and  doubly  ionized 
particles,  the  dimensionless  quantity,  7  has  been  calculated  by  taking  a  local  "Saha 
equilibrium”.  The  somewhat  involved  calculation  finally  leads  to  an  expression  which 
depends  on  Te  and  p  again.  The  expression  for  7  can  now  be  well  approximated  by 


7  ~  1.25  ■  10"5  1 


0.22/n 


•  t  +  0.85 


(A. 60) 


The  relative  error  of  thes  formula  amounts  to  about  ±5%  within  the  range  for  fr  and 
p  given  by  14000/f  <  ft  <  40000 A"  and  l(PPa  <  p  <  3  ■  105Pa,  respectively.  The 
maxium  electron  temperature,  Te,  within  the  core  of  an  argon  arc  follows  therefor  from 
the  eqs.(A.48)  and  (A.60)  by 

l 

rc 


Tc  = 


129  ■  Jo 


1  -  1.89  ■  10-3 {l  -  0.221n  (*)}  '  ^  *  (£) ' 


(A.61) 


where  Jq  is  slightly  varying  function  depending  on  n  and  £(rc)  and  given  according  to 
Table  (A.l). 

In  addition  to  the  temperature  maximum  Tc  within  the  arc  channel  of  the  constrictor 
one  is  interested  in  the  average  electron  temperature  defined  by 

f'  =  4  Ttrdr  (A. 62) 

r‘  Jt= 0 

Replacing  Te  in  the  integrand  by  eq.(A.42),  one  obtains 


=  TeJi  (A. 63) 
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where  the  coefficient  J\  is  a  function  of  the  current  density  profile  exponent,  n  and  the 
dimensionless  quantity,  £(rc)  «  (  «  £,  which  accounts  for  radiation  -  and  reaction  losses. 
In  Table  (A.'2)  this  coefficient  is  presented  for  various  values  of  n  and 


? 


S: 
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•  O 


Hydrogen 


0.01  bor 


Figure  A.3:  Plot  of  7  as  Funktion  of  te  and  p  as  Parameter  of  Hydrogen  as  Propellant 


l/r  [A/m] 

Figure  A.4:  Maximum  Electrontemperature  within  an  Hydrogen  Arc  Core  as  Function 
ot  the  Quotient  -j-  (arc  current  divided  by  the  core  radius)  for  differnt  pressure  p,n  = 
4,£=0. 
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